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'*^The  reduced  tolerance  Imaging  concept  Is  to  use  Imaging  system  hardware  of  reduced  complex¬ 
ity  to  make  phase-error-degraded  measurements  (or  to  lose  phase  Information  altogether)  anc 
then  reconstruct  diffraction- 11ml ted  Imagery  In  a  post- detection  processing  stage  using  a 
phase  retrieval  algorithm.  In  the  first  year  of  a  two-year  effort  several  advances  were 
made  toward  this  end.  An  estimation  theoretic  (Cramer-Rao)  lo»/er  bound  on  the  error  of 
estimating  a  coherent  Image  from  far-fleld  (Fourier)  Intensity  (squared  modulus)  measure¬ 
ments  was  derived  for  the  case  of  Gaussian  detector  noise.  Uniqueness  of  reconstruction 
from  Fourier  modulus  assuming  a  priori  known  support  was  proven  for  a  particular  class  of 
objects  --  sampled  objects  whose  support  area  outside  of  which  It  Is  zero)  has  a  con¬ 
vex  hull  with  no  parallel  sides.  A  closed-form  recursive  reconstruction  algorithm  was  de¬ 
veloped  for  reconstructing  such  objects  via  their  autocorrelation  functions.  Simulations 
showed  the  closed-form  solution  to  be  sensitive  to  noise  compared  with  iterative  Fourier 
transform  algorithms.  The  uniqueness  proof  was  useful  for  predicting  support  constraints 
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19.  for  which  iterative  reconstruction  is  facilitated.  Several  potential  constraints 
for  use  in  reconstruction  algorithms  were  examined  briefly,  but  support  and  non- 
negativity  are  the  only  two  constraints  that  have  been  extensively  exploited. 
Convergence  problems  when  the  support  constraint.  Imposed  on  the  world  by  active 
Illumination,  has  tapered  edges  were  ameliorated  by  a  modification  to  the  Itera¬ 
tive  transform  algorithm  using  an  "expanding  mask."  Alternative  reconstruction 
algorithms  were  studied.  Including  various  gradient  search  algorithms  (for  which 
analytic  expressions  for  the  gradient  of  the  error  metric  were  derived)  and  a 
modelling  approach,  but  they  have  not  yet  been  developed  to  the  point  where  they 
outperform  the  Iterative  transform  algorithm.  Laboratory  experiments  have  been 
planned,  starting  with  an  active  laser  Illumination  of  the  target  with  a  known 
Illumination  pattern  and  Fourier  Intensity  measurements*  Laboratory  experimental 
set  up  was  begun. 
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INTRODUCTION  AND  OVERVIEW 


1.1  BACKGROUND 

In  many  Imaging  scenarios  that  require  fine  resolution  at  long 
ranges,  phase  errors  limit  the  achievable  resolution  and  prevent 
diffraction-limited  Imaging.  The  phase  errors  may  arise  from  a  variety 
of  sources.  Including  atmospheric  turbulence,  misaligned  or  aberrated 
optics,  motion  compensation  errors,  local  oscillator  errors,  and 
waveform  generator  errors.  The  conventional  approach  for  obtaining 
diffraction- 11ml ted  Imagery  Is  to  build  Increasingly  more  complex  sensor 
hardware  having  tight  tolerances  on  Its  various  components  to  achieve 
the  desired  phase  stability. 

An  alternative  approach  Is  to  build  hardware  having  reduced 
tolerances  on  Its  phase  stability,  and  correct  for  the  phase  errors  by 
employing  a  phase  retrieval  algorithm  In  a  post-processing  step.  In 
some  Instances  a  sensor  can  be  used  that  Is  capable  of  measuring 
Intensity  only  and  does  not  measure  the  phase.  Then  a  phase  retrieval 
algorithm  Is  used  to  retrieve  the  lost  phase.  This  Is  what  we  refer  to 
as  Reduced  Tolerance  Imaging  (RTI).  Using  this  approach  one  can 
potentially  achieve  diffraction- 11ml ted  Imagery  using  a  sensor  system 
that  Is  less  complex,  cheaper,  lighter  weight  and  less  bulky. 

In  order  for  a  phase  retrieval  algorithm  to  work.  It  Is  necessary 
to  have  some  form  of  a  priori  Information  about,  or  constraints  on,  the 
Image.  Examples  of  such  constraints  that  have  been  useful  In  the  past 
are  nonnegativity  (applicable  to  Incoherent  Imaging)  and  knowledge  of 
the  object's  support  (knowing  Its  width  or  shape,  which  Is  available  for 
objects  on  dark  backgrounds  or  If  one  controls  the  pattern  of  radiation 
that  Illuminates  the  object). 
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Scvtral  liiportant  Issuts  mst  b«  addresstd  to  makt  th«  RTI  concept 
feasible.  Constraints  nust  be  found  that  art  poMerful  enough  to  ensure 
that  the  retrieved  phase  and  the  reconstructed  leage  are  uniquely 
related  to  the  eeasured  data.  The  relationship  between  the 
reconstructed  leage  and  the  Masured  data  eust  be  robust  enough  that  It 
Is  not  overly  sensitive  to  noise  or  other  laperfectlons  In  the  data  or 
constraints.  Reconstruction  algorithms  eust  be  found  that  converge 
reliably  to  a  solution  with  a  reasonable  amount  of  computation  and  In 
the  presence  of  realistic  amounts  of  noise. 

This  report  describes  the  results  of  the  first  year  of  a  two-year 
program  to  develop  the  Reduced  Tolerance  Imaging  concept. 

1.2  OVERVIEW  OF  ACCOMPLISHMENTS  TO  DATE 

In  this  section  the  principal  results  of  the  first  year  of  the  RTI 
program  w11'.  be  briefly  summarized.  They  are  reported  In  detail  In  the 
sections  and  appendices  that  follow. 

One  would  like  to  know  how  well  one  could  ever  hope  to  reconstruct 
an  Image  from  given  data  and  constraints.  Then  one  would  know  whether 
current  reconstruction  algorithms  are  good  enough  or  further  development 
Is  needed.  One  would  also  be  able  to  evaluate  and  compare  various 
measurement  schemes  without  having  to  develop  reconstruction  algorithms 
for  each.  This  can  be  done  using  estimation- theoretic  lower  bounds  on 
the  reconstruction  errors.  The  Cramer-Rao  lower  bound  was  derived  for 
the  case  of  far-fleld  Intensity  measurements  with  additive  Gaussian 
noise.  The  lower  bound  was  computed  and  compared  with  actual  errors 
experienced  In  Imagery  reconstructed  from  simulated  date.  These  results 
demonstrate  the  usefulness  of  estimation  theory  for  predicting  system 
performance.  Section  2  and  Appendix  0  describe  these  results. 

For  discrete*  or  sampled,  objects  of  a  certain  type  a  closed-form 
recursive  reconstruction  algorithm  has  been  developed.  It  reconstructs 
an  Image  from  the  autocorrelation  function  which  Is  gotten  by  Inverse 
Fourier  transforming  the  measured  Fourier  Intensity  data.  Although  the 
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closcd-for*  reconstruction  elgorltNi  has  questionable  usefulness  because 
It  Is  sensitive  to  noise.  It  has  provided  valuable  Insights  Into  the 
reconstruction  problea.  It  constitutes  a  uniqueness  proof  for  the  class 
of  objects  for  which  It  Is  applicable  and  suggests  lllualnatlon  pattern 
shapes  that  are  advantageous.  These  results  are  described  In  Section  3 
and  Appendices  A.  B  and  C. 

Since  leage  reconstruction  with  degraded  Fourier  phase  or  no 
Fourier  phase  requires  a  priori  constraints  on  the  object.  It  Is 
leperative  that  object  constraints  that  are  sufficiently  powerful  and 
robust  be  found.  The  vast  Majority  of  the  work  to  date  has  concentrated 
on  two  constraints:  support,  or  shape  (which  occurs  naturally  for 
iMgIng  satellites  and  aay  be  forced  by  an  llluelnatlon  pattern)  and 
nonnegativity  (valid  for  eost  passive  Incoherent  Imaging  probleas). 
Issues  relating  to  these  and  other  potential  constraints  are  discussed 
In  Section  4. 

Mhen  a  support  constraint  Is  leposed  by  using  an  active 
llliMlnatlon  pattern  at  the  target  to  achieve  the  desired  known  shape, 
the  principal  problea  Is  diffraction  effects  at  the  edges  of  the 
llluMlnatlon  pattern.  This  Makes  the  llluelnatlon  pattern  fall  off 
slowly  and  saoothly.  I.e..  Is  tapered,  rather  than  falling  off  abruptly 
as  would  be  preferred.  It  has  been  found  that  reconstruction  Is  Much 
easier  when  there  Is  little  or  no  tapering  of  the  llluelnatlon  pattern. 
Previous  versions  of  the  Iterative  reconstruction  algorlthe  were 
unsuccessful  In  reconstructing  coup  1ex> valued  Iwages  when  large  aeounts 
of  taper  was  present.  leproved  versions  of  the  algorlthe.  eeploylng  an 
"expanding  eask.”  were  developed,  and  this  resulted  In  a  greatly 
iMproved  result.  It  consists  of  using  an  unrealistically  seall  support 
constraint  for  early  Iterations,  which  forces  the  energy  of  the  Iwage  to 
be  better  centered  within  the  true  support  constraint,  and  using 
progressively  larger  support  constraints  for  later  Iterations.  Section 
5  describes  the  effects  of  different  types  of  llluMlnatlon  patterns, 
describes  the  Iwproved  algorlthe  eeploylng  the  expanding  eask.  and  shows 
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txptrlM*^U1  reconstruction  results. 

The  Iteretlvc  elgorlthi  described  In  Section  5  Is  one  of  several 
possible  approaches  to  solving  the  phase  retrieval  problee.  laproved 
ilgorlthas  are  sought  which  are  faster  and  aore  robust.  One  faally  of 
alternative  algorlthan  are  the  gradient  search  algorlthatt.  They  consist 
of  defining  a  aerlt  function,  coaputlng  the  gradient  of  the  aerlt 
function  as  a  function  of  a  paraaeter  space,  and  searching  In  the 
paraaeter  space  for  a  alnlaua  of  the  aerlt  function  In  the  direction  of 
the  negative  of  the  gradient  (the  global  alnlaua  of  the  aerlt  function 
defines  the  solution,  the  reconstructed  laage).  Merit  functions  that 
were  exaalned  Include  the  aaount  by  which  the  aodulus  of  the  Fourier 
transfora  of  an  object  estlaate  differs  froa  the  aeasured  Fourier 
aodulus  data  and  the  aaount  by  which  an  output  laage  violates  the 
object'doaain  constraints.  Paraaeter  spaces  that  were  Investigated 
Include  the  space  of  object  estlaates  and  the  space  of  Fourier  phase 
cstlaates.  Closed* fora  expressions  for  the  gradients  were  derived,  and 
the  entire  gradients  can  be  efficiently  coaputed  using  a  saall  nuaber  of 
fast  Fourier  transforas.  Gradient  search  algorlthas  were  tested  on  the 
coaputer  with  alxed  results  to  da^e,  but  they  show  proalse  and  will  be 
developed  further.  These  results  are  described  In  Section  6  and 
Appendices  E.  F  and  G. 

Another  approach  to  solving  the  phase  retrieval  problea  Is  a 
aodeling  approach.  The  coaplex  Fourier  transfora  or  pieces  of  It  are 
aodeled  by  a  paraaeterized  function.  The  aeasured  Fourier  aodulus  Is 
1east*squares  fitted  to  the  aodulus  of  the  aodel  to  deteralne  the 
unknoM  paraaeters.  Then  the  paraaeters  are  Inserted  Into  the  coaplex 
aodel  which  Is  evaluated  to  detenalne  the  phase.  Atteapts  to  aake  the 
aodeling  approach  work  were  unsuccessful.  It  Is  likely  that  the  aodels 
used  were  not  appropriate  to  the  coaplex  Fourier  transforas  of  Interest. 
Better  aodels  would  be  needed  before  further  work  along  these  lines 
should  be  pursued.  This  work  Is  discussed  In  Section  7. 
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Thf  vast  Mjorlty  of  tht  phasa  ratritval  work  prior  to  tha  currant 
effort  ravolvad  around  analysis  and  coaputar  slaulatlons.  Since  tha 
coaputar  slaulatlons  lapllcitly  assuaa  a  discrete  aodal  for  tha  object, 
there  Is  a  danger  that  the  real,  continuous  world  alght  behave 
differently.  For  this  and  other  reasons  It  Is  very  laportant  to 
deaonstrate  feasibility  on  real  data  collected  In  the  laboratory  that 
allows  us  to  Include  the  laportant  rea1*wor1d  effects  on  the  data.  At 
least  two  experlaents  will  be  perforaed:  an  active,  coherent  experlaent 
and  a  passive  Incoherent  experlaent.  The  active  coherent  experlaent  Is 
being  set  up  In  the  laboratory.  It  Includes  the  lllualnatlon  of  the 
target  with  a  laser  beaa  pattern  having  the  desired  lllualnatlon  shape 
and  controlled  aaounts  of  edge  tapering.  A  lens  foras  the  far-fleld 
(Fourier  transfora)  at  a  detector  plane.  A  second  channel  Including 
laaging  optics  will  be  used  to  fora  a  "ground  truth*  laage.  Section  8 
describes  the  active  coherent  experlaent  being  set  up  and  aentlons  plans 
for  the  passive  Incoherent  experlaent. 
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2 

INFORMATION  THEORETIC  LOMER  BOUND  FOR  PHASE  RETRIEVAL 
2.1  INTRODUCTION 

In  phase  retrieval  problems.  It  Is  desired  to  estimate  the  phase  of 
the  Fourier  transform  of  an  object  given  measurements  of  the  magnitude 
(I.e.,  the  modulus  or  the  square  root  of  the  Intensity)  of  the  Fourier 
transform.  This  Is  equivalent  to  estimating  the  object  Itself  because 
of  the  Fourier  transform  relationship.  Several  Iterative  Fourier 
transform  algorithms  have  had  great  success  In  making  such  object 
estimates  from  Fourier  magnitude  data  and  object  constraint  Information 
[2.1,  2.2].  However,  other  than  through  empirical  results  [2.3],  It  has 
not  been  known  how  the  error  In  the  object  estimate  depends  on 
measurement  noise,  constraint  Information,  and  other  parameters 
describing  the  problem. 

Results  In  estimation  theory  Include  a  number  of  methods  whereby 
lower  bounds  on  the  mean-squared  error  of  the  object  estimate  may  be 
calculated.  These  methods  use  knowledge  of  the  measurement  procedure, 
the  statistics  of  the  object,  and  the  statistics  of  the  noise  process  to 
compute  an  error  lower  bound.  An  Important  feature  Is  that  these 
methods  do  not  require  specification  of  the  algorithm  used  to  compute 
the  object  estimate  from  the  measured  data.  The  lower  bound,  then.  Is 
independent  of  the  algorithm  and  therefore  Indicative  of  the  best 
possible  estimation  performance  given  the  chosen  measurements  and  the 
underlying  statistics. 

The  Cramer-Rao  lower  bound  Is  the  most  often  used  lower  bound 
because  It  Is  usually  the  least  difficult  to  compute.  It  has  been  used 
In  a  large. number  of  single  and  multiple  parameter  and  time-varying 
waveform  estimation  problems  with  great  success  [2.4].  Algorithms  exist 
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which  produce  estimates  that  achieve  the  Cramer-Rao  bound  In  problems  In 
which  the  measurements  are  linearly  related  to  the  parameters  to  be 
estimated,  the  noise  Is  additive,  and  the  statistics  arc  Gaussian.  In 
nonlinear  problems  (of  which  phase  retrieval  will  be  an  example),  the 
lower  bound  can  usually  be  achieved  only  at  high  signal-to-nolse  ratios 
[2.4,  2.5];  nonetheless,  the  lower  bound  Is  generally  regarded  as  an 
Important  first  step  In  evaluating  and  designing  measurement  procedures 
and  parameter  estimation  algorithms  for  these  problems.  The  application 
of  lower  bounds  to  two-dimensional  signal  recovery  problems  described 
here  Is  a  recent  development,  and  It  Is  shown  that  It  Is  again  a  useful 
tool.  Appendix  D  gives  further  background  material  on  Cramer-Rao  lower 
bounds. 

2.2  PHASE  RETRIEVAL  PROBLEN  DEFINITION 

From  the  many  combinations  of  possible  phase  retrieval  problems  and 
underlying  assumptions,  the  following  specific  example  Is  chosen.  It  It 
desired  to  estimate  a  two-dimensional,  complex-valued  object  f.  from 

m 

real-valued  measurements  where  m  ■  (m^,  m^);  m^,  m^  *0,  1,...,  M-1 
and  p  "  (pj,  p^);  p^,  p^  "  0,  1,...,  2M-1.  The  measurements  are  related 
to  the  object  by 


and 


'pH  I  (2-2) 

m 

where  1^  Is  the  magnitude- squared  (Intensity)  of  the  discrete  Fourier 
transform  of  f,  c  Is  a  proportionality  constant,  N^  Is  additive  noise, 
<m,  p>  ■  m^p^  *  m^p^p  ind  summation  over  m  Implies  the  double  summation 
over  m^  and  m2.  Object  constraint  information  Is  essential  for 
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cstiMtIng  the  object.  The  weighting  array  w.  Is  explicitly  Included  In 
Eq.  (2-2)  to  allow  arbitrary  support  constraints  to  be  placed  on  the 
object.  For  an  object  of  M  by  M  resolution  elewents.  Nyquist  sampling 
requires  a  measurement  array  of  size  2M  by  2N  because  the 
magnitude- squared  has  twice  the  bandwidth  of  the  conq)1ex-va1ued  Fourier 
transform.  It  will  be  convenient  later  to  consider  w»  f,  S,  I.  and  N  as 
vectors.  The  phase  retrieval  problem  Is  to  estimate  the  object  f  given 
the  set  of  measurements  S  and  knowledge  of  the  constraint  that  the 
product  w^f^  Is  zero  wherever  w^  Is  known  to  be  zero. 


This  mathematical  statement  can  represent  a  number  of  applications 
In  which  phase  retrieval  problems  arise.  For  example*  consider  the 
measurement  geometry  shown  In  Fig.  2-1.  A  known*  complex-valued* 
monochromatic  wavefront  w(x*y)  Illuminates  an  unknown*  complex-valued 
object  f(x*y).  Alternatively*  for  the  wavefront  sensing  problem*  an 
unknown  monochromatic  wavefront  may  pass  through  a  known  aperture  having 
known  complex  transmittance  w(x*y).  The  Intensity  I(u*v)  In  a 
measurement  plane  located  a  distance  R  from  the  object  plane  Is: 


I(u*  v) 


(XR)^ 


I ][«(«.  y)f(x,  y)  .«p  dx  dy 


2 

(2-3) 


where  x  Is  the  wavelength  and  It  Is  assumed  that  R  Is  sufficiently  great 
that  the  Fraunhofer  approximation  can  be  made.  A  discrete  set  of 
measurements  S  Is  made  with 


Sp  ■  ,,T  J  I(u*  v)  du  dv  +  Np  (2-4) 
AA 

where  n  Is  the  detector  efficiency*  T  Is  the  detector  Integration  time* 
^A  Is  the  area  of  a  detector  element*  N.  Is  the  detector  noise*  and  the 
subscript  p  ■  (p^*  P2)  Indexes  over  the  measurement  plane.  A  phase 
retrieval  method  (e.g.*  an  Iterative  Fourier  transfona  algorithm)  would 
be  applied  to  the  measurement  set  S  using  the  object  constraint  provided 
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R - m 


Figure  2*1.  (U)  MeMurcmcm  Gcomciry  for  Phaic  Retrieval  ProMem 
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by  the  IHunlnatlon  pattern  w  to  give  an  estiMte  of  a  sanpled  version 
of  the  object  f.  Conversion  of  Eqs.  (2*3)  and  (2-4)  Into  discrete  form 
gives,  for  this  application,  a  value  for  the  constant  c  In  £q.  (2-1)  of 
nT^(Aa/XR)  where  Aa  Is  the  square  area  of  an  object  sample. 

The  complex- valued  object  f  can  be  written  In  terms  of  Its  real  and 
Imaginary  parts. 


f  ■  f*'*  -f 


Equation  (2-2)  then  becoaws 


I 


P 


2.3  CRAMER-RAO  LOWER  BOUND 


(2-5) 


(2-6) 


It  can  be  proven  that  the  variance  of  any  unbiased  estimate  of  a 
component  of  a  random  vector  Is  greater  than  or  equal  to  the 
corresponding  diagonal  element  of  the  Inverse  of  what  Is  called  the 
Fisher  Information  matrix.  The  value  of  the  diagonal  element  Is  the 
Cramer-Rao  lower  bound.  The  elements  of  the  Fisher  Information  matrix 
depend  in  turn  upon  the  second  partial  derivatives  of  the  Joint 
probability  distribution  of  the  measurement  vector  and  the  vector  to  be 
estimated.  This  result  Is  proven  primarily  by  the  use  of  the  Schwarz 
Inequality. 

Application  of  the  Cramer-Rao  method  for  determining  lower  bounds 
on  estimation  errors  to  a  specific  problem  must  therefore  begin  with  a 
determination  of  the  statistics  of  the  parameters  to  be  estimated  and  of 
the  noise  [2.4,  2.6].  In  this  analysis,  it  Is  assumed  that  f]!,  fj[,  and 
Np  are  each  statistically  Independent,  zero  mean,  Gaussian  random 
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variables  Mith  variances  Of/2*  end  respectively.  Note  that 

this  Implies  that  the  variance  of  the  comp  lex- valued  Is  Of. 

By  the  definition  of  conditional  probability* 

p(S.  f)  -  p(S|f)p(f)  (2-7) 

where  p(S.f)  Is  the  Joint  probability  density  of  S  and  f,  p(S|f)  Is  the 
conditional  probability  density  of  S  given  f»  and  p(f)  Is  the 
probability  density  of  f.  (Recall  that  f  and  S  are  vectors.)  The 
assumption  of  Gaussian  statistics  gives 


p(f)  "TT-t-  ®*p 

J.  J-tro^ 


m 


■ig' 


and,  using  Eqs.  (2-1)  and  (2-6)  which  Imply  that  p(S|f) 


P  n 


2o 


N 


(2-8) 

p(N  ■  S  -  cl), 

(2-9) 


The  Cramer-Rao  method  continues  by  defining  the  Fisher  Information 
matrix  J  In  terms  of  the  probability  density  functions.  For  the  present 
problem,  where  It  Is  desired  to  estimate  the  statistically  Independent 
real  and  Imaginary  parts  of  f,  a  workable  notation  Is  to  partition  J 
Into  four  submatrices: 


J 


[r  : 

• 

.j""  : 

(2-10) 


n 


2  2  2 

J  Is  of  dimension  2M  by  2N  (representing  the  M  Independent  plus 
2  i  I”  2 

the  H  Independent  fj)  and  each  of  the  submatrices  Is  of  dimension  M  by 

A  In 

The  elements  of  the  submatrices  are  defined  by.  for  example  [2.4, 

2.6], 


where  E[*]  denotes  expectation  taken  over  both  f  and  N  and  the  partial 
derivative  holds  S  constant.  The  other  submatrices  are  defined  by 
appropriate  substitution  of  the  superscripts  r  and  1.  It  Is  assumed  that 
these  and  any  other  required  derivatives  exist.  This  assumption  Is  valid 
for  the  phase  retrieval  problem. 

The  CrameroRao  method  concludes  by  determining  the  Inverse  of 
the  Fisher  Information  matrix  J.  The  diagonal  elements  of  are  the 
desired  lower  bounds  on  the  mean-squared  error  of  the  object  estimate  f. 

From  the  convention  used  to  define  J,  the  upper  left  diagonal  elements 

- 1  »•  •1-1 
of  J  refer  to  f„  and  the  lower  right  elements  to  f„.  If  J  Is 

M  fH 

similarly  partitioned  Into  four  submatrices: 


rK’’’’ :  K''n 


(2-12) 


then  the  Cramer-Rao  lower  bound  el_  on  the  mean-squared  error, 

A  A  Um 

E[|f  -  f_r].  In  the  estimate  f  of  object  component  f  .  Is  the  sum  of 
the  diagonal  elements  for  f|||  and  f|j^: 

■  C  *  C  (2-'3) 
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This  Is  the  quantity  which  the  following  analysis  seeks.  Strictly,  the 
lower  bound  Is  only  for  unbiased  estimates  of  f.  It  Is  beyond  the  scope 
of  this  work  to  determine  whether  particular  phase  retrieval  algorithms 
give  unbiased  estimates. 

2.4  LOWER  BOUND  FOR  PHASE  RETRIEVAL 


Substituting  Eqs.  (2-8)  and  (2-9)  Into  Eq.  (2-11),  differentiating, 
and  discarding  a  term  with  zero  expected  value  gives  [2.4] 


(2-14) 


Is  the  Kronecker  delta  function.  Similar  results  hold  for  the 

r1  1r 

Other  submatrices  of  J  except  that  J  and  J  have  no  term.  It  is 
Important  to  note  that  this  result  holds  for  any  function  I  of  the 
parameter  f.  It  does  not  assume  that  the  measurements  are  of  the 
Fourier  magnitude-squared. 


Equation  (2-6)  can  now  be  used  to  compute  the  first  term  on  the 
right  hand  side  of  Eq.  (2-14).  First, 


-iTr<J  -  m,  p> 


IT 


+  c.c. 


Then, 


(2-15) 


I 


m  n 


j  k 


exp 


-lTT<j  k  -  m  -  n,  p> 
M 
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“S-rZI  «"  -  ">  ^  n,  P>] 

j  k  *•  ■' 


+  c.c. 's 


(2-16) 


Taking  the  expected  value  gives 


(M7) 

P  J 

The  sumnatlon  over  k  Is  ellMlnated  since  the  and  fj[  are  Independent. 
The  first  and  third  terms  In  Eq.  (2>16)  are  also  eliminated  because  the 
and  fi  have  equal  variances.  Finally,  the  suimnatlon  over  p  gives 

m  nl 


because 


Z 


(2-18) 


(2-19) 


14 


Equation  (2-18)  Is  a  general  expression  for  one  of  the  submatrlces 
of  the  Fisher  Information  matrix  J  given  the  assumptions  above.  Similar 
computations  show  that  and  »  0.  In  this  case,  then, 
J  Is  diagonal  and  can  be  analytically  Inverted  to  obtain  This  Is, 
of  course,  a  result  of  the  discrete  Fourier  transform  nature  of  Eg. 
(2-2).  Other  phase  retrieval  problems  may  lead  to  nondlagonal  J 
matrices  which  may  be  difficult  or  Impractical  to  Invert  analytically. 


Using  Eqs.  (2-13)  and  (2-18),  the  lower 

mean-squared  error  In  the  estimate  of  f_  Is: 

ni 


2 

^f 


- T - 


1  + 


J 


bound 


on  the 


(2-20) 


It  Is,  as  stated  earlier.  Independent  of  the  phase  I'etrleval  algorithm 
used  to  estimate  f. 


The  notation  of  Eq.  (2-20)  can  be  simplified  by  defining  a 
signal-to-nolse  ratio: 


where,  by  Eq.  (2-6), 


SNR 


r 


ECelp]  ■  co|  J] 

J 

Equation  (2-20)  then  becomes 

2 

2  ^  _ ^f 

4  SNR  |w^  1^  * 

,  ^  'ml 
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(2-21) 


(2-22) 


(2-23) 


As  would  be  expected,  the  lower  bound  on  the  estiMte  reduces  to  the  a 
priori  variance  Of  either  Is  not  llliMlnated  •  0)  c.*  the  SNR 
Is  zero.  The  lower  bound  also  approaches  zero  as  the  SNR  ^;>proaches 
Infinity. 

For  the  case  In  which  the  magnitudes  of  the  support  constraint  w 
are  either  zero  or  one,  Eq.  (2-20)  predicts  that.  If  the  support 
constraint  Includes  a  smaller  part  of  the  M  by  M  object  array  (and 
thereforeZjwj  decreases),  then  the  error  lower  bound  Increases.  This 
Is  due  to'^he  loss  of  signal  as  can  be  seen  from  Eq.  (2-22).  On  the 
other  hand.  If  the  SNR  Is  held  constant,  then  Eq.  (2-23)  predicts  that 
the  error  bound  decreases.  This  Is  equivalent  to  sampling  at  greater 
than  the  Nyquist  rate  In  the  measurement  array  In  the  Fourier  domain. 
The  well-known  error  decrease  Is  known  as  compression  gain. 

It  Is  known  that  current  Iterative  phase  retrieval  algorithms  are 
more  successful  In  converging  to  a  solution  for  some  object  support 
constraints  than  for  others  (e.g.,  for  a  triangularly-shaped  pattern 
Imposed  by  w,  the  algorithm  more  readily  finds  a  solution  than  for  a 
square  pattern)  [2.7].  By  a  solution  is  meant  an  object  estimate  that 
Is  as  close  to  agreeing  with  the  measured  data  and  the  a  priori 
constraints  as  possible.  In  some  cases,  an  algorithm  stagnates  and 
produces  an  output  In  poor  agreement  with  the  data  and  constraints;  such 
an  output  should  not  be  considered  an  object  estimate.  If  there  Is  more 
than  one  solution  that  closely  agrees  with  the  data  and  constraints,  the 
algorithm  may  find  a  solution  that  Is  different  from  the  true  object. 
There  Is  a  tendency  for  Iterative  transform  algorithms  to  find  solutions 
more  readily  for  cases  guaranteed  to  have  unique  solutions  (e.g., 
objects  with  triangular  support  constraints).  However,  when  the 
solution  Is  unique.  It  Is  also  known  that.  If  a  solution  Is  found  (I.e., 
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the  algorithm  does  not  stagnate  In  poor  agreement  with  the  data  and 
constraints),  then  the  mean-squared  error  Is  Independent  of  the  shape  of 
the  object  support  constraint.  From  either  Eq.  (2-20)  or  Eq.  (2-23),  It 
can  be  seen  that,  for  a  given  value  ofZlwJ^,  the  lower  bound  e^ 
depends  only  on  |W||||  and  not  on  the  two-dimensional  distribution  of  w 
(the  support  constraint).  The  Cramer-Rao  lower  bound  Is  apparently  a 
measure  of  the  error  of  algorithms  which  have  found  a  reasonably  good 
estimate  and  Is  Insensitive  to  lack  of  uniqueness  or  to 
algorithm-dependent  problems  such  as  stagnation.  The  Insensitivity  to 
uniqueness  Is  further  demonstrated  by  an  example  shown  In  Appendix  D. 


2.5  CONCLUSION  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


In  this  Investigation  of  the  application  of  estimation  theoretic 
lower  bounds  to  phase  retrieval  and  image  reconstruction  problems,  the 
Cramer-Rao  lower  bound  on  the  mean-squared  error  In  the  object  estimate 
from  Fourier  magnitude-squared  measurements,  given  additive  noise, 
Gaussian  statistics,  and  Nyquist  sampling,  was  found.  The  lower  bound 
approaches  the  appropriate  values  In  the  limits  of  high  and  low  SNR,  but 
does  not  depend  on  the  object  support  constraint.  Further  research 
should  Investigate  other  measurement  models  (e.g.,  Fourier  magnitude 
measu'ements) ,  object  domain  constraints  (e.g.,  nonnegativity), 
statistical  assumptions  (e.g.,  Poisson  noise),  and/or  other  Information 
theoretic  lower  bounds  to  extend  and  refine  the  bounds  and  to  attempt  to 
show  a  dependence  on  a  priori  knowledge  such  as  support  constraints. 
Computer  simulations  and  laboratory  experiments  should  also  be  performed 
to  allow  comparison  of  the  lower  bound  to  the  error  achieved  by  current 
phase  retrieval  algorithms. 
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3 

UNIQUE  CLOSED  FORM  RECONSTRUCTION  ALGORITHM 
3.1  INTRODUCTION 

Since  the  object's  autocorrelation  function  can  be  computed  from 
the  modulus  of  its  Fourier  transform,  reconstructing  the  object  from 
its  autocorrelation  is  equivalent  to  reconstructing  it  from  the  modu¬ 
lus  of  its  Fourier  transform.  In  an  earlier  effort,  it  was  shown 
that  a  unique  closed-form  algorithm  for  reconstructing  an  object  from 
its  autocorrelation,  which  operated  in  a  recursive  fashion,  was  pos¬ 
sible  for  two  very  special  kinds  of  objects:  those  fitting  within  a 
rectangle  with  an  additional  point  off  one  corner  of  the  rectangle 
and  those  fitting  within  a  triangle  having  nonzero  corners.  This 
earlier  result  has  been  vastly  generalized  to  include  objects  having 
supports  whose  convex  hulls  have  no  parallel  sides,  a  very  large 
class  of  objects.  This  generalized  algorithm,  which  includes  a 
uniqueness  proof,  is  described  in  Section  3.2  and  Appendicies  A,  B 
and  C. 

Experimental  reconstruction  results  obtained  using  the  algorithm 
are  shown  in  Section  3.3.  Although  the  present  form  of  the  algorithm 
is  very  sensitive  to  noise,  limiting  its  practical  use,  it  has  proven 
to  be  very  valuable  in  that  it  suggests  useful  illumination  pattern 
(support)  constraints,  as  is  demonstrated  in  Section  4.1.  Another 
problen  with  this  reconstruction  algorithm  is  that  it  explicitly 
assumes  a  sampled  object,  i.e.  one  consisting  of  an  array  of  delta 
functions,  and  it  cannot  in  its  present  form  be  employed  for  real- 
world  continuous  objects.  One  possible  way  around  this  proolem  is  to 
use  the  quasi -samp ling  method  suggested  in  Section  3.4. 
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3.2  PHASE  RETRIEVAL  FOR  DISCRETE  FUNCTIONS  WITH  SUPPORT  CONSTRAINTS 


3.2.1.  INTRODUCTION 

The  reconstruction  of  object  functions  having  non-redundant 
spacings  was  discussed  In  [3.1].  Hayes  and  Quatlerl  [3.2]  showed 
that  the  boundaries  of  triangular  objects  can  be  reconstructed  by 
making  use  of  certain  spacings  In  the  object  which  are  non-redundant. 
In  another  direction,  Bruck  and  Sodin  [3.3]  showed  that  the  unique¬ 
ness  of  phase  retrieval  Is  equivalent  to  the  Irreduciblllty  of  a 
polynomial  In  two  variables  which  Is  closely  related  to  the  Fourier 
transform  (z-transform)  of  the  object  function.  Fiddy,  Brames  and 
Dainty  [3.4]  used  Elsenstlen's  Irreduciblllty  criterion  to  prove 
uniqueness  for  object  functions  satisfying  certain  support  con¬ 
straints  and  showed  that  Flenup's  Input-output  Iterative  Fourier 
transform  algorithm  [3. 5-3. 7]  converged  faster  to  a  better  recon¬ 
struction  when  these  constraints  were  satisfied.  Flenup  [3.8]  pre¬ 
sented  a  closed-form  algorithm  for  reconstructing  such  object  func¬ 
tions  from  their  autocorrelation  functions.  He  also  presen.-ed  a 
similar  closed-form  reconstruction  algorithm  for  objects  sat:sfying  a 
triangular  support  constraint  and  thereby  showed  that  such  objects 
are  uniquely  defined  by  their  autocorrelation  functions  among  all 
object  functions  satisfying  the  same  support  constraint. 

A  generalization  of  Flenup's  results  to  a  wider  class  of  support 
constraints  Is  presented  here.  Also,  an  algorithm  for  generating 
closed-form  reconstruction  algorithms  Is  described.  Brames  [3.9] 
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recently  obtained  a  result  similar  to  the  uniqueness  theorem  in 
Section  3.2.3. 

3.2.2.  MASKS 

2  2 
Let  A  denote  the  Euclidean  plane  and  let  Z  denote  the  points  in 

2  2 

with  integer  coordinates.  A  finite  subset  of  Z  is  a  mask  if  it 

contains  at  least  three  non-collinear  points  and  its  coiivex  hull  in 
^  (the  smallest  convex  set  containing  it)  has  no  parallel  sides. 

Let  M  be  a  mask  and  let  [M]  denote  its  convex  hull  inA  .  Then  [M] 
is  a  convex  polygon  (including  its  interior).  See  Figure  3-1.  A 
vertex  v  of  [M]  is  opposite  a  side  s  of  [M]  if  the  line  through  v  and 
parallel  to  s  contains  no  points  of  [N]  other  than  v  (see  Figure  3-2). 
A  vertex  of  [M]  is  a  reference  point  of  M  if  it  is  opposite  some  side 
of  [M]  (see  Figure  3-3).  The  set  of  all  reference  points  of  M  will  be 
denoted  by  R(M). 

3.2.3.  UNIQUENESS  THEOREM 

2 

Let  f  be  a  complex-valued  function  on  Z  .  The  support  of  a 
2 

function  on  Z  is  the  set  of  points  at  which  the  function  is  non¬ 
zero.  Let  *S(f)  denote  the  support  of  f.  If  »S(f)  is  a  finite  set, 

2 

the  autocorrelation  function  of  f  is  defined  for  x  c  Z  by 

r(x)  .  2^  f(y)f*(y  -  x)  (3-1 ) 

y«z^ 

where  the  *  denotes  complex  conjugation.  Let  f^  be  another 
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FIGURE  3-1.  IS  A  MASK.  M^  IS  NOT  A  MASK  SINCE  [M^]  HAS 
TWO  PARALLEL  SIDES. 
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FIGURE  3-3.  THE  CIRCLED  VERTICES  OF  [M]  ARE  THE  REFERENCE 
POINTS  OF  THE  MASK  H. 
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2 

complex-valued  function  on  Z  with  finite  support  *S(f^)  and 
autocorrelation  function  r^ . 

We  have  the  following  uniqueness  theorem. 

Theorem;  If  M  is  a  mask,  R(M)  9  j(f )  C  J(f  ^ )  C  M  and  r  ■  r^  then 

there  exists  a  complex  number  a  of  modulus  1  such  that  f^  -  af. 

The  proof  is  in  Appendix  A. 

3.2.4.  RECONSTRUCTION  ALGORITHMS 

In  this  section  closed-form  algorithms  for  reconstructing  a  func¬ 
tion  from  its  autocorrelation  function  will  be  described. 

Let  S  be  the  number  of  vertices  of  [M].  Let  Vq,  .  .  v^_^  be 

an  ordering  of  the  vertices  going  around  [M]  in  the  counter-clockwise 
direction  and  let  p^,  .  .  .,  pj_^  be  a  similar  ordering  of  the 
reference  points  of  M.  By  Lemma  A-2  In  Appendix  A,  R(M)  contains  an 
odd  number  of  points  so  that  T  is  odd.  Let  K  •  (T  -  l)/2  and  let 

Qn  ■  P(nK)niod  7  ^  ■  0»  •  •  •»  T  -  1.  Since  K  and  T  are 

relatively  prime,  the  are  distinct  and  hence  run  through  all 
the  points  of  R(M)  (see  Figure  3-4).  By  Lemma  A-4  in  Appendix  A, 

*^0  *^(n+l)mod  T  unique  separation  in  M.  That  is,  if  x, 
y  .  mnd  X  -  y  .  ^  th«  x  .  T 

y  ■  "n- 

Let  N  be  the  number  of  points  in  M.  A  reconstruction  algorithm 
for  the  mask  M  is  an  ordered  pair,  (q,  m),  where  q  -  (qg,  .  .  .,  q^^ 
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FIGURE  3-4.  THE  NUMBERINGS  OF  THE  VERTICES  AND 
A  MASK.  HERE  S  >  8.  T  -  5,  AND  K 


vj  *  P3  *  04 


2. 


POINTS  OF 


is  an  ordering  of  the  points  in  M  and  m  >  (myt  •  •  •»  a 

sequence  of  integers  satisfying  the  following  conditions.  The  points 

Qq,  •  •  •»  <lj  1  described  above.  For  n  -  T,  .  .  .,  N  -  1, 

the  integers  satisfy  the  conditions  0  £  £  T  -  1,  and 

M  n  (M  ♦  )  C  Iqq,  .  .  q^l  and 

n 

Mn(M  -  q^  ♦  q^  )  9  |  q^,  .  .  .,  q^_^  |  .  In  the  next  section,  an 
n 

algorithm  for  generating  such  reconstruction  algorithms  will  be  described. 

In  order  to  justify  the  above  definition  of  reconstruction  algo¬ 
rithms  it  will  now  be  shown  how  they  can  be  used  to  reconstruct  a 
function  from  its  autocorrelation  function. 

Let  f  be  a  complex-valued  function  on  satisfying  R(K)  9  S(f)  9  M 
and  let  r  be  its  autocorrelation  function.  Now  let 

*  ■  ^(n^l)mod  T  ■  ‘^n  suppose  that  for  some  y  c  Z^,  f(y)f*(y  -  x)  ^  0. 
Then  y  c  <S(f)  9  M  and  y  -  x  e  »S(f)9M.  Also,  y  -  (y  -  x)  ■  x  ■ 

"(n^Dmod  T  -  5I"«  «(„.,)„d  T 

separation  in  M,  it  follows  that  y  •  j  and  y  -  x  -  q^. 

2 

Therefore  y  -  j  the  only  y  e  Z  for  which 

f(y)f*(y  -  x)  ^  0,  hence 

'•<'’(ri*1)mod  T  -  On)  ■  '("(n.ljwd  t’  '*<"„).  (3-2) 

and  since  R(M)  9  *S(f)t  )mod  T  ~  ^  follows  from 

(3-2)  that 


|f(qo)l^ 


n!?0  '  ‘’{2n^l)mod  T^ 

- 

naU 


(3-3) 
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Since  f  is  defined  by  r  only  up  to  multiplication  by  a  modulus  1 
complex  number,  we  may  require  that  f(qQ)  >  0.  Then  TCQq)  is  equal 
to  the  positive  square  root  of  the  right-hand  side  of  (3-3).  Now 
f(q^)  can  be  computed  for  n  -  1,  .  .  T  -  1  from  the  formula 
f{q„)  -  '‘(qp  -  (%_})•  shown  in  Appendix  B  that  if 

(q,  m)  is  a  reconstruction  algorithm  then  the  following  program  will 

2 

compute  f(qp)  for  n  -  T,  ...»  N  -  1.  Set  f(x)  -  0  for  x  e  Z 
and  X  ^  q^,  n  -  0,  ....  T  -  1,  and  set  n  -  T  -  1. 

Step  1:  If  n  -  N  -  1,  stop.  Otherwise  n  n  ♦  1. 


Step  2: 


f(q„) 


-  "m  ) 
n 


n-1 

-  2  f(q|()  f*(q|, 

k-O 


q  ♦  q  ) 
n 


/f*(q 


m_ 


Step  3:  Go  to  Step  1. 


3.2.5.  ALGORITHM  FOR  GENERATING  RECONSTRUCTION  ALGORITHMS 

It  will  be  assumed  that  we  are  given  a  sequence  of  all  vertices 
''O’  '  ‘  *’  '^S-1  where  M  is  a  mask  and  the  sequence  is  ordered 

in  the  counter-clockwise  direction  around  [M]. 

For  n  ■  0,  .  .  . ,  S  -  1 ,  let  s^  be  the  side  of  [M]  with  end- 

2 

points  v^  and  ^jn+ij^iod  S’  ^  linear  operator  on  31 

which  rotates  each  vector  inJJ^  90*  counter-clockwise. 

First,  the  reference  points  p^,  .  .  .,  p^  ^  imtst  be  found.  Note 
that  every  side  of  [M]  has  a  vertex  opposite  it  which  is  therefore  a 
reference  point.  Of  course,  several  sides  may  have  the  same  vertex 
opposite  them.  Let  w^  -  5  -  v^  for  n  -  0 . S  -  1. 
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A  vertex  is  opposite  a  side  s^  if  and  only  if 

^'^n^  for  k  -  0,.  .  S  -  1,  where  <,>  denotes 
the  inner  product  on  (see  Figure  3-5).  Thus,  by  taking  each  side 

in  the  order  Sq . .  all  the  reference  points  of  M  will  be 

found,  and  if  they  are  numbered,  in  the  order  in  which  they  are  found, 
Pq,  .  .  .,  Pj  then  the  ordering  will  be  in  the  counter-clockwise 
direction  around  [M]. 


As  mentioned  above  T  Is  odd.  Let  K  ■  (T  -  l)/2  and  - 

P(nK)niod  T*  ”  *  ^  Since  each  is  a  reference  point 

and  therefore  is  a  vertex  of  [M],  there  exists  an  Integer  k^  such  that 


0  £  k^  £  S  -  1  and  q„  -  v^  .  For  n  -  0,  .  .  .,  T  -  1,  define 


J'n  -  “"(VD-OX  S  -  t' 


(3-4) 


Then  by  Lenina  A-3  in  Appendix  A,  for  x  c  M,  x  4  q„  and  x  4  *1(0+1  )(nod  T’ 
yo>  <  <x,  y„>  <  <q(o+i)TOd  T»  J'n^*  equivalent  to 

saying  all  points  in  M  excluding  q^  and  q(o+i)inod  T  strictly  between 
lines  pendicular  to  y^  and  passing  through  q^  and  ^(n+ijioojj  y  See 
Figure  3-6.  (The  uniqueness  of  separation  of  q^^  and  q(o+i)nK)d  T 
mentioned  in  Section  3.2.4  follows  from  this  double  inequality.) 


NO.  let  .  q„  ♦  T  "X  l«t  .  .„/2  for 

n  ■  0,  .  .  .,  T  -  1.  Then  is  the  midpoint  of  the  line  segment 
joining  q^  and  q(n+ij,„Q(j  y  Let  0  -  M  ^  R(M)  (set  difference) 

2 

and  let  #  be  the  characteristic  function  of  0  as  a  subset  of  Z  . 

2 

This  is,  4  is  the  function  on  Z  which  is  1  on  D  and  0  outside  D. 
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figure  3-5. 


AN  ILLUSTRATION  OF  THE  VECTORS  AND  Uw^.  HERE  n  -  0 
AND  THE  ORIGIN  IH  §?  IS  DENOTED  BY  "0". 
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For  X  c  D  and  0  £  n  £  T  -  1 ,  let  -  <  x  -  >  •  Then 

'hj^(x)  I  / 1  ly^l  I  is  the  distance  from  x  to  the  line  through  8^  and 
perpendicular  to  y^,  txhere  Mypll  denotes  the  Euclidean  norm  of  y^ 
(see  Figure  3-7).  The  set  0  contains  N  -  T  points  and  we  define  T 
orderings  of  the  points  in  0,  0  •  I  q.  •  •  *.  |  . 

n  -  0 . T  -  1,  satisfying  >  I  ^  ^ 

k  -  0,  .  .  N  -  T  -2.  The  following  program  generates  sequences 

f  *  *t  ^Nl  ^  ^^T  ^  ^  1  * 

Set  naT-1  andk-Oand  enter  the  following  loop. 

Step  1:  If  n  -  N  -  1,  stop.  Otherwise  define 

b  •  min  jj:0<j£N-T-land  #(d|^  j)  ■  1  }  . 

Step  2:  If  -  cl,^  Ij)  •  1.  go  to  Step  7. 

Step  3:  n  •  n  ♦  1 . 

Step  4:  Define  •  d. 

n  K  ,0 

step  5:  If  h.  (q,)  >  0,  define  -  k.  Otherwise  define 
m^  .  (k  1  )mod  T. 

Step  6:  #(q^)  •  0. 

Step  7:  k  -  (k  ♦  1  )mod  T  and  go  to  Step  1. 

It  is  shown  in  Appendix  C  that  the  loop  is  not  infinite  and  if 
q  ■  (Pq,  ....  and  m  •  (niy,  .  .  .,  )  Then  (q,  m)  is  a 

reconstruction  algorithm. 
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3.2.6.  IMPLEMENTATION 


The  algorithms  presented  in  the  last  two  sections  can  be  imple¬ 
mented  with  two  computer  programs.  The  first  program  would  implement 
the  algorithm  in  Section  3.2.5.  Its  input  would  be  a  mask  and  its 
output  would  be  a  reconstruction  algorithm.  The  second  program  would 
implement  the  program  in  Section  3.2.4.  Its  input  would  consist  of  i 
reconstruction  algorithm  and  an  autocorrelation  function  and  its  out¬ 
put  would  be  the  object  function.  With  this  arrangement,  if  one 
wished  to  reconstruct  many  object  functions  using  the  same  mask,  the 
first  program  would  have  to  be  run  only  once. 

3.2.7.  CONCLUSIONS 

It  has  been  shown  that  if  a  function  is  zero  outside  a  given  mask 
and  is  non-zero  at  the  reference  points  of  the  mask,  then  it  is 
uniquely  determined  (up  to  multiplication  by  a  complex  number  with 
modulus  1)  by  its  autocorrelation  function  among  all  other  object 
functions  which  are  zero  outside  the  mask.  (A  mask  is  a  set  of 
points  in  the  discrete  lattice  whose  convex  hull  has  no  parallel 
sides.)  Moreover,  there  is  an  algorithm  for  generating 
reconstruction  algorithms  for  any  given  mask  which  in  turn  can  be 
used  to  reconstruct  object  functions  satisfying  the  above  mentioned 
conditions  from  their  autocorrelation  functions. 

This  theory  has  some  similarity  to  holography  [3.10,  3.11]. 
However,  here  several  (at  least  3)  reference  points  are  used  whereas 
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only  one  reference  point  Is  needed  in  the  holographic  situation.  On 
the  other  hand,  the  holographic  reference  point  must  be  isolated  from 
the  rest  of  the  object  whereas  no  such  isolation  of  the  reference 
points  is  required  here.  It  is  interesting  to  speculate  whether 
there  might  be  a  more  general  theory  of  which  this  theory  and  holo¬ 
graphy  would  both  be  special  cases. 


35 


3.3  EXPERIf€NTAL  CLOSED-FORM  RECONSTRUCTION  RESULTS 

Autocorrelation  data  was  computer-simulated,  including  the 
effects  of  noise,  and  images  were  reconstructed  using  the  closed-form 
reconstruction  algorithm  described  in  the  previous  section. 

I 

For  each  reconstruction  experiment  an  object,  f(x,y),  fitting 
within  a  triangular  support  was  Fourier  transformed: 

F(u,v)  -^f(x,y)] 

The  squared  modulus,  |f(u,v)|^,  of  the  Fourier  transform  was  com¬ 
puted,  and  it  was  scaled  in  intensity  so  that  the  total  integrated 
intensity  became  equal  to  a  given  number  of  photons, 

uv 

Then  each  intensity  sample  F(u,v)r  was  replaced  with  a  random 

12  ' 

drawn  from  a  Poisson  distribution  with  mean 
12  12 

and  variance  equal  to  F(u,v)  .  When  F(u,v)  is  a  large  num¬ 
ber  {>32)t  then  a  Gaussian  approximation  to  the  Poisson  distribution 
is  used.  This  Poisson  noise  process  simulates  the  effect  of  photon 
(shot)  noise  on  the  measured  Fourier  intensity  data.  The  normalized 
RMS  error  (NRMSE)  of  the  data  is  given  by 


Data  NRMSE  - 

A  noisy  autocorrelation  was  computed: 

r„{x,y)  j  ; 

and  an  image,  gp(x,y),  was  reconstructed  using  the  closed-form 
reconstruction  algorithm.  The  NRMSE  of  the  reconstructed  image  is 
given  by 
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Image  NRMSE 


E 

|ag^(x,y)  -  f(x,y)p 

|f(x.y)|^ 

xy 

where  a  Is  a  constant  chosen  to  minimize  the  error  metric,  which 
accounts  for  the  unknown  phase  constant  associated  with  f(x,y).  It 
can  be  shown  that  the  value  of  a  that  optimizes  the  Image  NRMSE  Is 


a 


f(x.y)g*(x.y) 

xy 


Examples  of  Images  of  objects  reconstructed  from  noisy  data,  for 
which  the  object  Is  a  uniform  triangle,  are  shown  In  Figures  3-6  to 
3-10  for  various  sizes  of  the  object.  Figure  3-11  -plots  the  image 
NRMSE  versus  the  data  NRMSE  for  the  images  shown  In  Figures  3-6  and 
3-9.  Several  interesting  effects  are  evidenced  from  these  recon¬ 
struction  examples.  First,  the  closed-form  algorithm  is  very  sensi¬ 
tive  to  noise.  A  fraction  of  a  percent  error  In  the  data  results  in 
several  percent  error  In  the  image.  Second,  increased  data  error 
results  in  increased  Image  error,  but  only  In  an  average  sense. 
Occasionally  the  Image  error  can  be  greater  when  the  data  error  is 
less,  because  for  a  given  amount  of  data  error  the  Image  error  that 
one  gets  is  highly  variable.  Depending  on  the  particular  realization 
of  the  noise  in  the  data,  the  three  corner  points  will  have  different 
amounts  of  error.  Small  differences  in  the  error  of  the  corner 
points  can  yield  large  differences  in  the  emr  of  the  image  since 
the  corner  points  are  used  over  and  over  again  and  the  error  from 
them  propagates  and  Is  magnified  as  the  recursive  steps  build  upon 
one  another.  This  also  gives  rise  to  a  third  effect:  the  error  for 
the  interior  points  of  the  reconstructed  image  is  much  worse  than  the 
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FIGURt  3-H.  8  .  8  OBJLC'  AND  IMAGLS  Rt CONSTRUCTED  BY  TFFE  CLOSED-FORM 

RECURSIVE  ALA'RrtF’-V  NuDbpr  of  photons  listed  are  in  the  Fourier  intensity 
doma 1 n. 
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FIGURE  3-9.  16  x  16  TRIANGULAR  OBJECT  AND  IMAGES  RECONSTRULTED  BY  TtIF 

CLOSED-FORM  RECURSIVE  ALGORITHM 


FIGURE  3-10.  32  x  32  TRIANGULAR  OBJECT  AND  IMAGES  RECONSTRUCTED  BY  THE 

CLOSED-FORM  RECURSIVE  ALGORITHM 


hJaPIMflLIZEO  F1M3E  IN  TME  ORTR 


FIGURE  3-11.  NRMS  ERROR  OF  THE  RECONSTRUCTED  IMAGE  VERSUS  NRMS  ERROR  OF 
THE  DATA.  Crosses  are  for  the  8x8  triangles  and  squares  for  the  16  x 
16  triangles. 
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error  of  the  edge  points  In  the  Image.  Fourth,  for  similar  reasons 
the  error  of  larger  images  is  far  worse  than  the  error  of  smaller 
images.  Fifth,  a  stripe  pattern  parallel  to  each  of  the  edges  tends 
to  occur.  This  happens  for  the  following  reason.  Suppose  that  one 
of  the  three  reconstructed  corner  points  is  brighter  than  it  should 
be.  Then  the  opposing  edge  of  the  image,  the  computation  of  which 
involves  division  by ttrc corner  point,  will  tend  to  be  too  dark.  Then 
the  next  inward  row  (or  column)  from  the  edge,  the  computation  of 
which  Involves  subtraction  of  terms  involving  the  edge,  will  tend  to 
be  too  bright,  etc. 

How  badly  this  striping  effect  affects  the  interpret ability  of  an 
image  was  tested  by  using  a  picture  of  an  airplane  as  the  object, 
imbedded  in  a  32  by  32  triangle.  The  results  of  reconstruction 
experiments  from  noisy  data  using  the  closed-form  reconstruction 
algorithm  for  the  object  are  shown  in  Figures  3-12,  3-13  and  3-14. 
As  seen  from  Figure  3-12,  the  image  can  still  be  discerned  through 
the  partially-obscuring  striped  pattern.  Therefore  the  intelligi¬ 
bility  of  the  image  may  be  understated  by  the  image  NRMSE. 
Figure  3-13  shows  the  image  NRMSE  versus  the  total  number  of  photons 
for  all  the  noise  values  tried  for  this  object,  and  Figure  3-14  shows 
the  same  information,  but  as  a  function  of  data  NRMSE. 
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FIGURE  3-12.  32  x  32  JET  OBJECT  AND  IMAGES  RECONSTRUCTED  BY  THE  CLOSED- 

FORM  RECURSIVE  ALGORITHM.  The  number  of  photons  in  the  Fourier  intensity 
data,  the  corresponding  data  NRMS  error  and  the  reconstructed  image 
NRMS  errors  are  as  follows: 


No.  of  Photons 

Data  NRMSE 

Imaqe  NRMSE 

(a) 

(b) 

q  -  Orlg 

10®  , 

Inal  Object  - 
0.0011 

0.03 

(c) 

0.0034 

0.10 

(d) 

6  X  IO7 

0.0043 

0.11 

(e) 

2  X  10' 

10' 

0.0073 

0.27 

(f) 

0.0107 

0.60 
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IMAGE  NRMSE 


INAGE  NRMSE 


1.0 

0.8 


0.4 


DATA  NRMSE 


FIGURE  3-14.  RECONSTRUC'n’D  JET  IMAGE  NRNS  ERROR  VERSUS  DATA  NRMS  ERROR 
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3.4  QUASI -SAMPLING  ILLUMINATION  PATTERN 

One  probleM  with  the  closed-fone  reconstruction  algorithm 
described  in  Section  3.2  is  its  reliance  on  a  sampled  object.  In 
this  section,  it  is  shotm  that  by  a  special  kind  of  illumination  one 
can  approximately  achieve  the  desired  sampling  effect  in  the  target 
area. 

If  one  illuminates  a  target  area  Kith  four  mutually  coherent 
point  sources  in  the  far  field  at  a  distance  R  given  by 

<(u  -  Uq.  *  *(u  -  Ug,  V  ♦  v^)  ♦  i(u  ♦  Uq,  V  -  Vq) 

♦  4(u  ♦  Uq  V  ♦  Vq) 

one  gets  a  field  at  the  target  area  given  by  the  sum  of  four  plane 
Kaves  (Fourier  transforms  of  the  oelta  functions): 


which  has  Intensity 

,.,™- ...» c-y) 

which  has  lines  of  zeros  along  x  .  xR{n  ♦  1/2)/(2Uq)  and  along 

y  .  xR(n  ♦  1/2)/(2Uq)^  for  n  -  0,  *1.  *2 . This  is 

illust*  ted  in  Figure  3-15. 
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FIGURE  3-15.  QUASI-SAMPLING  ILLUMINATION  PATTERN.  The  circles  are  the 
locus  of  zero?  ^ 
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Tm  lllualnatlon  pattern  would  be  of  llaUtd  extent  which  could  be 
Modeled  by  Multiplying  the  above  by  a  slowly  varying  weighting  function 
defining  the  f leld-of-vlew.  (XR/4)t(x,y ) ,  so  that  the  entire 
llluMlnatlon  pattern  Is 

.(x.y)  •  cot  cot 

Where  there  are  a  nueiber  of  cycles  of  the  cosines  over  the  extent  of 
t(x,y). 

What  Is  accomplished  by  this  Is  a  quasl-sampi Ing  of  the  object. 

It  may  be  possible  to  use  the  closed-form  recursive  reconstruction 
algorithm  to  reconstruct  an  object.  Illuminated  by  w(x,y)  above, 
from  Its  autocorrelation,  or  it  may  be  necessary  to  make 
modifications  to  reduce  the  errors  due  to  approximating  this  pattern 
by  a  true  sampling  pattern. 

Note  that  by  the  addition  of  more  plane  waves  It  is  possible  to 
qet  sharper  local  maxima  and  broader  stripes  of  low  intensity,  but 
at  the  expense  of  a  more  complicated  illumination  system  with 
phase-stability  problems  of  its  own. 

In  the  real  world,  the  four  mutually  coherent  illumination 
sources  may  have  an  unknown  relative  phasing  between  them.  If  the 
constant  relative  phases  of  the  four  sources  arm  ^  ^2,  ^3  and 
^4,  then  the  product  of  cosines  like  the  equation  above  occurs 

only  if  #1  -  #2  -  #3  ♦  <4  ■  0.  This  implies  a  stringent 
stability  requirement  on  the  illumination  system,  but  it  requires 

the  control  of  only  a  single  parameter  (one  piston  term)  rather  than 
the  control  of  the  phase  of  an  entire  large  aperture. 
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4 

CONSnUINT  INVESTI6ATI0N 


In  this  snctlofi  tiM  various  form  of  constraints  that  sight  ba  usaful 
for  phasa  ratrlaval  ara  dlscussad.  Sactlon  4.1  dascribas  rasults 
obtalnad  with  a  varlaty  of  support  (llluslnatlon  pattam)  constraints. 
Tha  vast  sajorlty  of  tha  affort  to  data  concantratad  on  daval oping 
Isaglng  concapts  basad  on  tha  support  constraint.  Sactlon  4.2  dascribas 
othar  constraints  that  sight  also  prova  usaful. 

4.1  EFFECT  OF  ILLUMINATION  PATTERN  SNAPE 

Tha  support,  S,  of  an  objact  Is  dafinad  as  tha  sat  of  points  for 
which  tha  objact  Is  nonzaro.  For  tha  casa  of  a  satalllta  Isagad  against 
tha  night  sky  or  a  ship  Isagad  on  cals  watar  with  a  SAR,  tha  support  of 
tha  objact  Is  basically  tha  f11  lad-in  outllna  of  tha  objact.  For  in 
airborna  or  spacaborna  sansor  looking  downward  at  a  ganaral  scana,  tha 
axtant  of  tha  objact  Is  basically  dafinad  by  tha  flald-of-vlaw  of  the 
sansor.  This  lattar  casa  does  not  raprasant  a  usaful  support 
constraint.  Howavar,  for  an  Isaglng  systas  asploying  active 
llluslnatlon,  tha  transsittad  baas  (tha  llluslnatlon  baas)  can  taka  on  a 
known  shape  at  tha  plana  of  tha  target,  and  It  can  ba  dasignad  to  occupy 
an  area  ssallar  than  tha  flald-of-vlaw  of  tha  racalvar.  Than  the 
affective  support  of  tha  objact  Is  tha  support  of  tha  llluslnatlon  baas 
pattam.  For  tha  casa  of  a  SAR,  It  Is  assusad  that  whan  no  phasa  Is 
available  tha  pulsa  'apatitlon  fraquancy  Is  at  least  twice  that 
ordinarily  required  by  lyquist  saspllng  whan  phasa  Inforsatlon  is 
available. 

Tha  two  sost  Isportant  properties  of  an  llluslnatlon  pattern  ara 
Its  shape  (elliptical,  rectangular,  polygonal,  ate.)  and  Its  taper  (how 
slowly  It  transitions  fros  tha  bright  part  of  tha  pattam  to  where  it  Is 
affectively  zero).  As  shown  In  tha  proposal  [4.1,  p.  2-29],  phase 
retrieval  algorithm  ara  such  sore  affective  for  sosa  shapes  (which  wa 
refer  to  as  strong  shapes)  than  for  others.  Furtharsora,  phase 
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rttritval  algorlthas  art  tort  affactivt  for  sharp  support  constraints, 
l.t.  whtn  thtrt  Is  lUtlt  or  no  taptr  to  tht  llliaalnaclon  pattarn  [4.1, 
p.  2*25].  Stctlon  S  of  this  rtport  dttalls  tht  rtsults  of  our 
Invtstigatlon  of  tht  tfftcts  of  taptrtd  llluailnatlon  patttrns  and 
algorltha  laprovtatnts  that  Mrt  aadt  for  tht  cast  of  largtr  atounts  of 
taptr.  In  what  lantdlattly  follows  wt  discuss  tht  tfftcts  of  tht  shape 
of  tht  support. 

In  tarly  phast  rttritval  work  thtrt  was  not  an  awartntss  that  tht 
support  of  tht  ohjtct  playtd  an  Important  rolt  In  tht  succtss  of  phase 
retrieval.  Early  successful  reconstruction  rtsults  were  for  space 
objects  whose  supports  wart  naturally  non-centrosyantrlc  [4.2].  Other 
groups  atttaipttd  phast  retrieval  for  unnatural  objects  —  scents  bounded 
by  squares  ••  and  were  unsuccessful.  Fiddy,  Brants  and  Dainty  [4.3] 
found  that  the  Iterative  Fourier  transfoni  algorithm,  although  it  worked 
poorly  for  a  rectangular  support,  worked  well  for  a  support  consisting 
of  a  rectangle  plus  an  extra  point  Just  off  one  comer  of  tht  rectangle. 
This  latter  support  has  tht  special  property  that  any  sampled  function 
defined  on  that  support,  which  Is  nonzero  at  the  extra  point  and  at  one 
opposite  corner,  has  a  Fourier  transform  that  Is  a  nonfactorable 
polynomial  according  to  Elsensteln's  Irreduciblllty  theorem.  This 
Implies  that  the  phase  retrieval  problem  In  unconditionally  unique  for 
objects  of  this  type.  In  retrospect,  from  those  results  we  can  make  the 
crucial  connection  between  three  different  aspects  of  the  phase 
retrieval  problem:  the  support  of  the  object,  the  uniqueness  of  phase 
retrieval,  and  the  success  of  the  Iterative  Fourier  transform 
algorithms. 

The  trends  connecting  those  three  elements,  which  we  have  continued 
to  confirm,  are  the  following.  First,  the  support  of  the  object 
determines  whether  ambiguities  are  possible.  Second,  objects  for  which 
uniqueness  can  be  proven  are  easier  to  reconstruct  by  the  Iterative 
Fourier  transform  algorithm  than  are  other  objects.  The  first  trend  Is 
amply  demonstrated  In  Section  3.2  which  shows  that  sampled  objects 
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having  known,  convax  hulls  with  no  parallal  sidts  art  uniqua.  Tht 
stcond  trtnd  Is  shown  by  tht  rtconstnictlon  rtsults  [4.1.  pp.  2-24  and 
2-28]  In  which  objtcts  having  known  triangular  support  (which  art  uniqut 
—  stt  Stctlon  3.2)  and  objtcts  having  known  supports  with  saparattd 
parts  (which  tvtn  In  ont  dlatnslon  art  usually  uniqut  ••  stt  txaapits  in 
Stctlon  5)  art  taslly  rtconstructtd  whilt  objtcts  with  othtr  support 
constraints.  Ilka  that  of  a  singit  tlllpst  or  a  singit  rtctangit  art 
difficult  to  rtconstruct. 

Tht  clostd'fona  raconstructlon  algorltha  dtscribtd  In  Stctlon  3.2 
■ay  not  bt  practical  for  ust  on  rta1*wor1d  data,  sinct  It  rtquirts  tht 
objtcts  to  bt  BOdtlltd  discrtttly  (at  a  grid  of  dolta  functions  or 
saapltd  points)  and  It  Is  vary  stnsitivt  to  nolst.  particularly  If  tht 
vtrttx  points  art  dia  (stt  Stctlon  3.3).  Ntvtrthtitss.  It  dots 
constltuta  a  uniqutntss  proof  for  tht  typts  of  objtcts  to  which  it  can 
In  thtory  bt  appiltd:  objtcts  whost  support  has  a  convax  hull  with  no 
parallal  sidts.  This  Itads  us  to  considtr  lllualnatlon  patttms  of  this 
typt.  Figurt  4-1  shows  an  txaaplt  of  a  raconstructlon  txptrlatnt  using 
an  lllualnatl on-pa tttm  shapt  suggtsttd  by  tht  uniqutntss  proof.  On  the 
Itft  Is  tht  Modulus  of  a  conpitx-valutd  SEASAT  SAR  lingt  Mil  tipi  ltd  by  i 
binary  pattam  (rtprtstnting  tht  IlluMlnatlon  patttm)  In  tht  shapt  of  a 
ptntagon.  In  tht  can  tar  Is  tht  Modulus  of  Its  Fourltr  transfora  (tht 
Fourltr  phast  was  discardtd).  Tht  Itaratlva  Fourltr  transfoni  algorlthM 
was  ustd  to  rtconstruct  an  iMagt.  tht  Modulus  of  which  Is  shown  on  tht 
right.  froM  tht  Fourltr  Modulus  using  tht  known  support  patttrn.  The 
result  shown  Is  after  several  hundred  Iterations  (It  had  not  coMpletely 
converged  yet),  and  It  strongly  rtstoblts  the  original  object,  although 
not  perfectly.  Given  the  difficulty  In  reconstructing  coMpltx- valued 
Inagts  with  contiguous  supports  (with  the  exception  of  triangular 
support)  [4.1.  pp.  2-24  to  2-30].  the  success  of  this  kind  of  support 
constraint  would  have  been  difficult  to  anticipate  wart  It  not  for  the 
uniqueness  proof. 
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FIGURE  4-1.  RECONSTRUCTION  EXPERIMENT  USING  PENTAGON-SHAPED 
ILLUMINATION  PATTERN,  Left  to  right:  modulus  of  complex-valued 
Illuminated  SEASAT  SAR  scene;  modulus  of  signal  history  (Fourier 
transform);  modulus  of  Image  reconstructed  from  modulus  of  signal 
history  and  pentagon-shaped  support  constraint  using  the  Iterative 
Fourier  transform  algorithm. 
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In  sumry*  to  dato  tha  support  constraints  found  to  bo  aost  ustful 
for  phasa  ratrlaval  ara  (a)  supports  having  two  or  aora  wall-saparatad 
parts  and  (b)  supports  having  convax  hulls  with  no  parallal  sidas.  Tha 
farthar  tha  sidas  ara  froa  balny  parallal »  tha  battar. 

4.2  OTHER  CONSTRAINTS 

Any  InfouMtlon  that  Is  In  a  doaain  othar  than  tha  doaain  of  tha 
aaasurad  data  has  tha  potantlal  for  baing  a  usaful  constraint  for  phasa 
ratrlaval.  Constraints  In  tha  doaain  of  tha  Baasurad  data  usually  just 
Halt  tha  avallabla  data  and  tand  not  to  ba  usaful  for  phasa  ratrlaval. 
Candidata  constraints  for  tha  phasa  ratrlaval  probloa  ara  Mstad  In 
Tabla  4.1. 


Tabla  4-1 

CANDIDATE  CONSTRAINTS 

Support  (lllualnatlon  pattam) 

Nonnagativlty 

Polarization 

Transaittad  wavafora 

Point  scattarar  In  scana 

Othar  scana  charactarl sties 


Of  thasa,  tha  support  constraint  racalvad  tha  aost  attantlon.  and  It  Is 
discussad  In  Sactlon  4.1.  Tha  othar  constraints  ara  dascribad  In  what 
follows. 

Nonnagativlty 

Tha  nonnagativlty  constraint  has  baan  vary  usaful  In  pravlous  phasa 
ratrlaval  afforts  [4.2]  and  also  In  othar  laaga  raconstructlon  problaas. 
such  as  toaographic  raconstructlon  froa  Incoaplata  projactlons  and 
constralnad  deconvolution.  Unllka  tha  support  constraint,  nonnagativlty 


55 


■ust  txist  naturally  —  m  do  not  know  how  to  lapost  It  artificially. 
It  naturally  occurs  In  aost  passive,  noncoheront  laaglng  scenarios.  The 
brightness  distribution  of  an  Incoherently-llluBlnated  reflecting  objact 
or  a  se1f-«n1ss1ve  object  Is  characterized  by  a  real,  nonnegatlva 
function  (power  or  photons  per  unit  area).  This  can  be  true  for 
actively  lllualnated  objects  as  well,  as  long  as  the  lllualnatlon  Is 
sufficiently  Incoherent.  An  exception  in  which  this  constraint  nay  not 
bfi  valid  Is  for  passive  doppler  laeglng  as  encountered  In  the  PACE 
p 'ogran  [4.4].  for  which  the  aperture  function  Is  band-pass. 
f/Onnegat1v1ty  Is  not  as  useful  for  band-pass  systeas  since  the  Inpulse 
response  has  very  large  negative  (or  coap lex- valued)  sldelobes  which 
convolve  the  Inage.  destroying  Its  nonnegativity.  For  the  cases  In 
which  nonnegativity  naturally  does  occur.  It  should  be  relied  on  heavily 
as  a  phase  retrieval  constraint. 

Polarization 

Certain  kinds  of  reflecting  objects  have  distinctly  different 
reflectivities  for  the  two  different  receive  polarities  (I.e.  for  sane 
polarity  as  transnit  or  for  opposite  polarity).  As  an  exanple,  comer 
reflectors  reflect  either  very  strongly  or  very  weakly  depending  on  the 
polarization.  Unfortunately,  froai  a  single  collection  It  Is  not 
iHMdlately  obvious  how  to  utilize  this  Infomatlon.  On  the  other  hand, 
If  two  collections  are  nade  sinultaneously,  one  for  each  polarization, 
then  there  Is  Increased  possibility  of  using  polarization 
advantageously.  One  such  possibility  would  be  to  use  the  difference 
between  two  degraded  Inages  (with  neasured  Fourier  phase  In  the  presence 
of  phase  errors)  to  Identify  point-llke  reflectors  Then  the  point-llke 
reflectors  could  be  used  In  the  proninent-point  processing  described 
below.  Other  exanples  of  using  polarization  nay  be  also  be  possible. 
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Transaitttd  Wivtfoni  Typt 

Early  on  In  tht  prograa  It  ms  thought  that  parhaps  transalttlng  a 
pul  St  with  Missing  frtgutncy  bands  sight  bt  ustful  Insofar  as  It  would 
constltuto  a  support-likt  constraint  In  tht  signal  history.  Upon 
furthtr  txaailnatlon.  It  apptars  that  such  nlsslng  frtqutnclts  would 
prlaarlly  rtsult  In  a  loss  of  data  rathtr  than  constituting  a  ustful 
constraint. 

A  point  worth  Mking  rtlating  to  transaltttd  wavtfront  Is  that  tht 
ust  of  phast  rttritval  ttchnlquts  aay  facllltatt  tht  ust  of 
nonconvtntlonal  Mvtforas.  As  tht  transaiitttd  Mvtfont  dtparts  froa  tht 
standard  stt  (t.g.  tht  chirp  wavtfora),  tht  availability  of  hardMrt 
that  can  fora  tht  dtsirtd  wavtfora  In  a  phast-stabit  aanntr  aay  bt 
qutstlonablt.  By  rtducing  tht  toltranct  on  tht  phast  stability  of  a 
wavtfora  gtntrator  it  may  be  possible  to  achieve  waveforms  that  would 
otherwise  bt  very  difficult  to  product.  Tht  reduced  tolerance 
iMgIng/phast  rttritval  ttchnlquts  nay  provide  the  ntans  for  rtducing 
the  phast  stability  of  tht  wavefora  gtntrator  while  Mintalning  the 
dtsirtd  resolution. 

Point  Scatttrtrs  In  Scent 

Presently,  point-1 Ikt  scatttrtrs  (pronintnt  points)  In  tht  target 
area  art  used  for  correcting  saull  anounts  of  phast  errors  In  SAR  signal 
histories  [4.5,  4.6].  Pronintnt  point  processing  can  also  be  of  great 
utility  for  the  cast  of  severe  phast  errors  or  when  no  phast  Infonution 
at  all  Is  neasured.  One  particular  scenario  for  phast  correction  In  the 
presence  of  large  ont-dlntnslonal  phast  errors  has  already  been 
denonstrattd  [4.6].  For  the  cast  of  notion  conpensatlon  errors  In  SAR, 
one  has  a  ont-dlntnslonal  (azlnuth)  phast  error.  This  occurs 
particularly  for  the  cast  of  Inverse  SAR,  for  txanpit  tht  radar  is 
ground-based  and  the  (noncooptrativt)  target  flits  by  with  a  poorly 
known  flight  path  and  rotation.  If  there  exists  a  doninant  proninent 
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point  scatttrtr  In  a  givtn  coaprtssod  rang#  cell,  th«n  It  can  be  used  to 
calculate  the  azleuth  phase  error  (taking  Its  phase  to  be  the  phase 
error).  The  phase  errors  In  all  range  bins  can  be  corrected  by 
subtracting  that  phase. 

Other  Scene  Characteristics 

The  constraints  Mentioned  above  are  coMMon  to  large  classes  of 
iMagery.  Also,  there  May  often  be  additional  constraints  that  exist  Ip 
specific  Instances.  For  exaeple.  If  the  scene  has  been  iHaged  by 
another  sensor  systeM  or  by  a  siMlIar  sensor  at  an  earlier  tlee,  then 
these  additional  leages  May  contain  Infonaatlon  that  can  be  counted  on 
to  appear  In  the  present  leage  and  therefore  can  be  used  as  an  a  priori 
constraint.  Exaeples  Include  the  known  existence  of  peraanent  cultural 
targets  or  of  no*retum  areas  such  as  lakes  or  saooth  surfaces. 
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RECONSTRUCTION  OF  OBJECTS  WITH  TAPERED  ILLUMINATION 


5.1  statement  of  PROBLEM 

It  Is  well  known  that  knowledge  of  the  support  of  an  object  can  be 
a  powerful  source  of  inforaatlon  in  image- reconstruction  problems.  By 
support  we  mean  a  compact  region  outside  of  which  the  object  Is  known  to 
be  zero,  and  we  denote  the  set  of  points  that  make  up  the  support  by  the 
symbol  S.  In  particular,  considerable  success  >^as  been  realized  In 
reconstructing  an  object  from  Its  Fourier  modulus  and  a  known  support 
[5. 1,5. 2].  In  the  reduced-tolerance  Imaging  program  an  effort  Is  being 
made  to  exploit  this  ability. 

Consider  an  active  sensor  system  that  illuminates  a  target  area  so 
that  the  illumination  Is  confined  to  a  predetermined  region.  Let  h(x,y) 
be  the  complex  reflectivity  of  the  target: 

h(x,y)  •  h(x,y)|  (5-1) 

Let  w(x,y)  be  the  complex  Illumination  function: 

w(x,y)  -  iw(x,y)t  (5-2) 

He  define  the  effective  object  as  the  product  of  the  complex 
reflectivity  of  the  target  and  the  Illumination  function: 

f(x,y)  «  w(x,y)  h(x,y) 

•  If(x,y)|  (5-3) 
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The  effective  object  Mill  now  have  a  support  corresponding  to  the  known 
extent  of  w(x,y).  The  Intensity  pattern  of  the  field  esanatli.g  fron  the 
llluMlnated  target  Is  Masured  In  the  far  field  which  aay  be  Interpreted 
as  the  squared  awdulus  of  the  Fourier  transfone  of  the  effective  object. 
Known  phase-retrieval  algorlthas  nay  then  be  eaiployed  to  reconstruct  the 
effective  object  froe  the  support  constraint  and  the  Measured  Fourier 
modulus. 

Notice  that  there  Is  some  freedom  In  the  choice  of  the  form  of  the 
Illumination  pattern.  For  example,  the  shape  of  the  outline  of  the 
pattern  could  be  specifically  selected  to  enhance  the  usefulness  of  the 
support  constraint.  It  Is  known  that  certain  symmetries  in  object 
support  can  create  stagnation  problems  In  phase- retrieval  algorithms. 
Consequently  the  outline  of  the  Illumination  pattern  should  have  an 
asymmetric  shape.  Furthermore,  there  Is  some  evidence  that  a  support 
consisting  of  disjoint  regions  can  be  an  advantage  In  phase  retrieval. 
Finally,  It  Is  useful  to  choose  an  Illumination  function  with  a  constant 
modulus  over  most  of  the  region  of  illumination  thus  facilitating  the 
Inversion  of  Eq,  (5-3): 


h(x,y) 


(x.y)eS 


i(0{x,y)-(t)  (x.y)) 

w 


(5-4) 


when  one  desires  the  complex  reflectivity  of  the  target  without  the 
Influence  of  the  illumination  pattern. 

Unfortunately,  the  modulus  of  the  Illumination  pattern  will  not  be 
binary  In  practice,  but  will  have  some  taper  associated  wit  at  the 
edges,  due  to  the  effects  of  diffraction  by  the  apertur-t^  of  the 
Illuminator.  The  contrast  between  an  Ideal  untapered  Illumination 
pattern  and  a  more  realistic  illumination  function  Is  Illustrated  In 


61 


T 


Figurt  5-1.  Intuitively  one  eight  expect,  end  experleentally  It  has 
been  shown  [5. 1,5. 2],  that  the  reconstruction  of  an  object  from  Us 
Fourier  mhIuIus  and  support  would  be  aort  challenging  for  an  object  with 
a  tapered  profile  than  for  one  with  a  sharp  profile.  It  was  the  purpose 
of  this  Inquiry  to  explore  this  Issue  via  coeputer  sleulatlon  and  to 
look  for  algorlthalc  aodlflcatlons  that  would  enhance  restoration  for 
this  case.  For  exaeple.  It  was  hoped  at  the  outset  that  any 
difficulties  Incurred  by  tapered  lllual nation  eight  be  offset  by  the 
support  being  disjoint. 

5.2  PRELIMINARY  SIMULATIONS 

Ue  began  by  exploring  the  effect  of  tapered  llluelnatlon  on  phase 
retrieval  by  eeans  of  coeputer  sleulatlon.  A  pair  of  disjoint  ellipses 
was  used  as  the  basic  shape  for  the  llluelnatlon  pattern.  The  untapered 
llluelnatlon  pattern  was  assigned  a  value  of  unity  within  the  ellipses 
and  zero  outside.  Taper  was  Introduced  by  convolving  the  binary 
ellipses  with  a  convolution  kernel.  The  norealized  kernels  used  in 
these  prellelnary  sleulatlons  are  shown  In  Figure  5-2.  Cross  sections 
of  the  edge  of  the  resulting  llluelnatlon  patterns  are  given  In  Figure 
5-3. 

As  aientloned  earlier.  It  was  speculated  that  disconnected  support 
eight  help  to  overcoee  any  problees  associated  with  tapered 
llluelnatlon.  For  this  reason  the  total  llluelnatlon  pattern  was  chosen 
to  be  two  disjoint  ellipses.  Sleulatlons  were  perforeed  for  objects 
with  differing  aeounts  of  llluelnatlon  taper  and  differing  amounts  of 
separation  between  ellipses  In  the  llluelnatlon  pattern.  A  given 
sleulatlon  was  perforeed  by  first  eultiplying  cojplex  SEASAT  SAR  Imagery 
by  the  given  llluelnatlon  pattern  to  create  an  effective  object. 
Because  It  Is  the  effective  object  that  we  try  to  recover  through 
phase-retrieval  techniques,  we  will  henceforth  refer  to  this  as  the  true 
object.  This  object  was  Fourier  transfomed  with  an  FFT  and  the  Fourier 
magnitude  was  retained.  The  known  region  of  support  In  the  object 
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FIGURE  5-2.  DISCRETE  CONVOLUTION  KERNELS  USED  TO  ADD  TAPER 
TO  BINARY  ILLUMINATION  PATTERN.  A.  Ccntcr-MighUd  ktmtl 
yields  taper  #1.  B.  Evenly-weighted  kernel  yields  taper  #2. 
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FIJURE  5-3.  CROSS  SECTIONS  OF  ILLUMINATION  OF  TAPER  USED  IN 
PRELIMINARV  SIMULATIONS 


dOMln  MS  supplied  by  herd  Halting  (thresholding)  the  lllueilnatlon 
pattern  with  a  very  smII  threshold  value.  A  standardized  sequence  of 
error-reduction  and  hybrid  Input-output  Iterations  [5.3]  were  then 
perforaed  to  reconstruct  the  object  froa  Its  Fourier  aodulus  and 
support.  Convergence  for  all  slaulatlons  ms  aonitored  by  calculating  a 
Fourier  doaain  noraallzed  error  aetric. 


E 


2 

F 


^  ('G(u.v)'  -  ,F(u.v)|)^ 

^ _ 

^  lF(u.v)l^ 

u.v 


(5-5) 


where  F(u.v)  Is  the  discrete  Fourier  transform  of  the  true  object  f(x,y) 
and  6(u*v)  Is  the  Fourier  transform  of  the  Image  estimate.  The 
convergence  Is  portrayed  In  Figure  5-4  for  six  kinds  of  Illumination- 
three  amounts  of  taper,  each  with  two  amounts  of  separation  between 
ellipses.  It  Is  Important  to  note  that  Figure  5-4  Is  a  log-log  plot  and 
therefore  the  behavior  of  the  algorithm  becomes  horizontally  compressed 
with  Increasing  number  of  Iterations.  Figures  5-5,  5-6,  and  5-7  give 
the  final  reconstructions  for  each  of  the  cases  tested.  These  results 
confirm  our  expectation  that  Increased  amounts  of  Illumination  taper 
make  the  reconstruction  process  more  difficult.  In  fact,  for  the  case 
with  the  largest  amount  of  taper  the  algorithm  convergence  appears  to 
have  stagnated.  This  Is  In  spite  of  the  fact  that  the  amount  of  taper 
Is  extremely  mild.  There  are  51  pixels  along  the  major  axis  of  the 
large  ellipse  and  only  two  pixels  of  taper  at  the  edge.  Thus 
convergence  appears  to  be  relatively  sensitive  to  Illumination  taper. 
It  Is  Important  to  realize  that  the  convergence  curves  shown  In  Figure 
5-4  correspond  to  a  specific  Initial  estimate  and  that  the  convergence 
behavior  could  vary  when  alternative  Initial  estimates  are  used. 


5.3  THE  SHRUNKEN-MASK  ALGORITHM 


In  order  to  explore  the  reasons  for  stagnation  we  created  a 
difference  Image  between  the  modulus  of  the  true  object  and  that  of  the 
restored  object  for  the  case  of  Intermediate  taper  (taper  #1).  This 
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FIGURE  5-4.  CONVERGENCE  BEHAVIOR  AS  A  FUNCTION  OF  ILLUMINATION  TAPER 
AND  SUPPORT  SEPARATION 
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FIGURE  5-5.  RECONSTRUCTIONS  OF  OBJECTS  WITH  UNTAPERED  ILLUMINATION 
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FIGURE  5-6.  RECONSTRUCTIONS  OF  OBJECTS  WITH  MILDLY  TAPERED 
ILLUMINATION.  (Taper  #1  In  Figure  5-3) 
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FIGURE  5-7.  RECONSTRUCTIONS  OF  OBJECTS  WITH  TAPERED  ILLUMINATION. 
(Taper  #2  In  Figure  5-3) 
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differenct  iMgt  Is  bipolar  and  a  bias  was  added  for  display  In  Figure 
5-8.  Notice  that  the  difference  Inage  Indicates  that  the  reconstruction 
Is  shifted  In  the  horizontal  direction  relative  to  the  true  object. 
This  suggests  that  the  algorltha  nay  be  stagnating  because  of  Its 
Inability  to  properly  register  the  reconstruction  relative  to  the 
support  constraint  when  tapered  llluni nation  Is  used. 

To  better  understand  this  conjectured  node  of  stagnation  consider 
an  object  with  tapered  llluninatlon  f(x»y).  We  define  a  binary  mask 
n(x.y)  that  Is  the  characteristic  function  of  the  known  support: 


m(x,y) 


1  ,  (x,y)eS 
0,  (x,y)eS' 


(5-6) 


where  S'  stands  for  the  conplenent  of  S.  An  Inage  g'(x»y)  outputted 
by  the  Iterative  Fourier  transform  algorithm  Is  the  Inverse  Fourier 
transform  of  a  Fourler-domain  estimate  having  modulus  equal  to  the  given 
Fourier  modulus  data  coupled  with  the  current  estimate  of  the  Fourier 
phase.  Suppose  the  output  image  Is  Just  a  shifted  version  of  the 
object: 


g'(x,y)  -  f(x  -  Xq,  y  -  y^) 


(5-7) 


A  shift  In  the  object  domain  Introduces  a  linear  phase  factor  in  the 
Fourier  domain  and  has  no  effect  on  the  Fourier  modulus.  This  output 
Image  will  clearly  satisfy  the  Fourier  modulus  constraint.  The  output 
''mage  has,  however,  been  shifted  relative  to  the  mask  so  that  the  object 
lain  support  constraint  has  been  violated.  In  other  words, 
.«jlt1ply1ng  by  the  mask  function  will  crop  an  edge  of  the  output  Image. 
We  use  a  normalized  error  metric  to  Indicate  the  degree  of  Inconsistency 
between  an  estimate  and  the  object  support  constraint: 


71 


FIGURE  5-8.  MODULUS  DIFFERENCE  BETWEEN  OBJECT  AND  RECONSTRUCTION. 
(Illumination  due  to  Taper  #?  in  Figure  5-3)  The  bipolar  difference 
image  has  been  biased  up  for  display. 
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(5-8) 


E 


2 

0 


|g'(x.y)m'(x,y)  ^ 

^  |g'(x.y)l^ 
x.y 


where  m'(x,y)  Is  the  cheracterlstic  function  of  S'.  If  the  shift  vector 
(Xo^yo)  Is  small  with  respect  to  the  Illumination  taper  the  object 
domain  error  metric  will  also  be  relatively  small.  This  is  because  only 
the  tapered  edges,  where  there  Is  little  energy,  will  be  cropped  and 
this  contributes  to  only  a  small  portion  of  the  total  object  energy. 

Though  the  cropped  output  Image  now  satisfies  the  support 
constraint.  Its  Four1er«transfonn  modulus  no  longer  exactly  equals 
|F(u,v)|  .  It  can  easily  be  shown  that  the  Fourler>doaMln  error  metric 
Is  also  small.  Thus  the  error  metric  penalty  is  small  in  either  domain 
when  shifting  a  tapered  object  by  a  small  amount.  An  algorithm  that 
chooses  successive  estimates  based  upon  these  error  metric  objective 
functions  will  be  Insensitive  to  small  shifts  and  would  easily  stagnate 
due  to  extremely  small  slopes  In  the  objective  function.  Such  an 
algorithm  would  be  ineffective  at  finding  the  proper  object 
registration.  Furthermore,  one  can  Imagine  that  with  the  right 
redistribution  of  the  cropped  object  energy  an  object  estimate  could 
correspond  to  a  local  minimum  In  the  objective  function. 

Although  the  mode  of  stagnation  Just  presented  is  conjecture,  it 
provides  the  motivation  for  the  "shrunken-mask"  algorithm.  The 
shrunken-mask  algorithm  Is  designed  to  find  the  proper  registration 
early  on  In  the  Iterative  reconstruction  thus  circumventing 
shift-related  stagnation  that  might  otherwise  appear. 

Consider  a  new  binary  mask  m^(x,y)  created  by  hardlimiting  the 
tapered  illumination  function  with  some  Intermediate  threshold  value: 
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m^(x,y)  = 


1  ,  (x.y)  such  that  w(x,y)|  >  t 
0,  (x.y)  such  that  w(x,y),  _<  t 


(5-9) 

where  t  1s  the  threshold  value.  0  t^  1.  Notice  that  «^(x,y)  will  be 
a  “shrunken"  version  of  the  full  mask  m(x,y)  defined  for  t  »  0.  Suppose 
that  we  employ  the  shrunken  mask  as  the  support  constraint.  If  we  crop 
the  true  object  with  the  shrunken  mask  this  will  yield  an  estimate  with 
a  modest  penalty  In  both  the  object  and  Fou*‘1er  domains,  so  long  as  the 
threshold  value  is  not  too  large.  Notice,  however,  that  a  shift  In  this 
cropped  estimate  will  yield  an  object-domain  penalty,  due  to  the 
shrunken  mask  and  the  artifically  created  discontinuous  object  et  es. 
that  IS  much  greater  than  the  penalty  that  would  be  due  to  tn*  normal 
support  constraint.  Tiius  one  would  expect  the  output  image  to  be 
centered  better  with  the  shrunken  mask. 

While  the  Fou  ler  modulus  and  the  shrunken-mask  support  construints 
are  inconsistent,  they  may  still  be  Jointly  enforced  in  an  Iterative 
reconstruction  algorithm  to  get  an  Intermediate  reconstruction.  We 
might  expect  this  Intermediate  result  to  display  gross  features  of  the 
true  Object  in  proper  registration.  Enlarging  the  mask  to  its  full  size 
(setting  t  *  0)  removes  the  constraint  Inconsistency  and  allows  for  a 
complete  reconstruction  that  hopefully  avoids  shift-related  local 
minima.  The  shrunken-mask  algorithm  Is  shown  schematically  in  Figure 
5-9. 

The  shrunken-mask  algorithm  was  first  tested  on  the  elliptical 
object  where  the  taper  (taper  #2  in  Figure  5-3)  induced  stagnation  in 
previous  trials.  The  convergence  characteristics  are  displayed  in 
Figure  5-10.  It  Is  clear  that  the  conventional  algorltnm  performed 
better  early  on  In  the  Iterative  sequence.  This  Is  reasonable  since  the 
support  constraint  Is  Initially  looser  and  easier  to  satisfy.  By 
contrast,  the  shrunken-mask  algorithm  error  metric  quickly  levels  off 
while  an  Intermediate  reconstruction  Is  being  produced  but  drops 
dramatically  when  the  full-size  mask  Is  Introduced. 
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FIGURE  5-9.  THE  SHRUNKEN-MASK  ALGORITHM 
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FIGURE  5-10.  CONVERGENCE  FOR  SHRUNKEN-KASK  ALGORITHM.  Taper  Is  Taper 
#2  In  Figure  5*3.  The  shrunken  nask  had  a  threshold  value  t  ■  .9. 
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5.4  THE  ENLARGING  MASK  ALGORITHM 

While  the  success  In  the  shrunken-iMSk  algorithm  Is  encouraging  the 
amount  of  Illumination  taper  for  which  It  worked  remains  extremely 
small.  A  much  more  substantial  taper  was  introduced  by  using  a  circular 
convolution  kernel  with  a  radius  of  4  pixels.  The  resultant 
illumination  pattern  is  shown  In  Figure  5-11.  When  the  shrunken-mask 
algorithm  was  applied  to  an  object  with  this  Illumination  the 
convergence  was  not  much  better  than  the  conventional  algorithm.  This 
was  true  for  a  variety  of  threshold  values  that  were  tested.  Apparently 
the  increased  taper  is  a  significant  obstacle  for  the  shrunken-mask 
algorithm. 

Recall  that  the  shrunken-mask  algorithm  Jumps  from  a  small  mask  to 
the  full  mask  In  a  single  step.  A  logical  generalization  of  the 
shrunken-mask  algorithm  uses  several  intermediate-size  masks  in  order  to 
make  a  more  gradual  transition  to  the  full  size  mask.  We  call  this  the 
"enlarglng-mask"  algorithm.  The  collection  of  masks  used  in  a  given 
application  Is  characterized  by  a  sequence  of  threshold  val«cs.  Tne 
convergence  curve  for  the  enlarging-mask  algorithm  when  applied  to  an 
object  with  this  increased  taper  is  shown  In  Figure  5-12.  The  scallop 
effect  exhibited  by  the  convergence  curve  is  due  to  the  successive 
application  of  increasingly  enlarged  masks.  The  enlarging-mask 
algorithm  clearly  out-performs  the  shrunken-mask  algorithm  and  the  final 
reconstruction  exhibits  very  gooj  agreement  with  the  data  and  support 
constraint. 

A  final  trial  was  performed  with  an  even  more  realistic 
Illumination  taoer  created  with  a  Gauss1an-1  Ike  convolution  kernel  with 
a  maximum  radius  of  6  pixels.  This  kernel,  K(r),  was  formed  by 
correlating  a  circle  function  with  a  radius  of  2  pixels  with  Its  own 
autocorrelation; 
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(5-10) 


where  the  double  star  Indicates  two-dlaenslonal  crosscorrelation.  Thl^ 
kernel  Is  a  close  approxlaatlon  to  a  two>dlMns1ona1  Gaussian  function. 
The  resultant  lllualnatlon  pattern  Is  shown  In  Figure  5*11.  Note  that 
this  lllualnatlon  has  a  saoother  taper  and  that  the  tails  extend  out 
further  at  very  low  levels.  The  convergence  curves  for  this  case  are 
shown  In  Figure  5-13.  Again  the  enlarglng-aask  algorltha  succeeds  at 
finding  a  reconstruction  that  Is  In  excellent  agreeaent  with  the  data 
and  support  constraint  whereas  the  conventional  algorltha  did  not.  This 
reconstruction  Is  visibly  Indistinguishable  froa  the  true  object.  The 
results  of  reconstructions  perforaad  with  and  without  the  enlarglng-aask 
algorltha  are  given  In  Figure  5-14  for  1  lliailnatlon  patterns  due  to  the 
circular  and  Gaussian  convolution  kernels. 

While  tapered  lllualnatlon  presents  significant  stagnation  probicas 
for  conventional  phase-retrieval  algorlthas.  these  exaaples  deaonstrate 
that  the  enlarglng-aask  algorltha  successfully  clrcuavents  these 
difficulties,  even  In  the  presence  of  large  aaounts  of  taper. 
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FIGURE  5-14.  RECONSTRUCTIONS  Wllh  AND  WITHOUT  THE  ENLARG I NG- MAS K 
ALGORITHM  (EMA),  The  1 1 1  uinl na 1 1  on  pattern  In  A-C  1s  shown  1n  Figure 
5-llB.  The  Illumination  pattern  In  0-F  Is  shown  In  Figure  5-llC. 
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6 

GRADIENT-SEARCH  METHODS  IN  PHASE  RETRIEVAL 


6.1  INTRODUCTION 

Researchers  have  explored  many  approaches  to  solving  the  phase 
retrieval  problem.  These  Include  direct  nethods  using  complex  zeros  in 
the  analytically  extended  Fourier  modulus  [6.1],  the  error-reduction 
algorithm  [6.2,  6.3],  input-output  algorithms  [6.3],  recursive 
algorithms  [6.4,  6.5],  and  gradient-search  algorithms  [6.3,  6.6,  6.7, 
6.8].  Of  these  approaches  the  Input-output  algorithms  or  more 
specifically  the  hybrid  Input-output  (HID)  algorithm  appears  to  be  the 
current  algorithm  of  choice  when  operating  on  2-d1mens1ona1  data.  The 
HIO  algorithm  has  consistently  outperformed  competing  algorithms  with 
respect  to  computational  burden  and  robustness  to  noise.  In  spite  of 
the  relative  success  of  the  input-output  algorithms  there  are  documented 
instances  In  which  such  an  algorithm  converges  extremely  slowly  or  even 
stagnates  in  its  convergence  [6.9]. 

In  this  report  we  are  interested  in  the  specific  phase- retneva I 
proDlem  for  which  the  Fourier  modulus  and  an  object  support  constraint 
are  knovwi.  We  resurrect  the  idea  of  employing  a  gradient-search  method 
In  the  hopes  of  developing  an  algorithm  that  will  compete  well  with  or 
complement  the  Input-output  approach.  Gradient-search  approaches 
require  the  determination  of  an  objective  function  that  indicates  the 
degree  of  consistency  with  the  data  and  the  constraints.  This  choice  is 
pivotal  In  designing  a  specific  gradient-search  algorithm.  We  propose 
here  three  distinct  objective  functions  and  explore  the  performance  of 
each  when  used  in  conjunction  with  standard  gradient-search  techniques. 
In  the  next  section  we  discuss  the  error-reduction  algorithm,  the  parent 
of  the  Input-output  algorithms  and  Indicate  how  It  can  be  Interpreted  as 
a  gradient-search  algorithm.  We  introduce  the  first  new  objective 
function,  called  the  summed  objective  function.  In  Section  6.3.  T;,e 
second  and  third  objective  functions  are  introduced  in  Sections  6.4  and 
6.5.  These  objective  functions  utilize  the  same  object-support  error 


83 


metric  but  differ  In  their  underlying  parameters.  We  conclude  in 
Section  6.6  with  projections  of  future  work. 

6.2  THE  ERROR-REDUCTION  ALGORITHM 


An  iterative  algorithm  that  has  enjoyed  much  success  in  phase 
retrieval  is  known  as  the  error- reduction  (ER)  algorithm,  which  may  be 
easily  understood  by  referring  to  Figure  6-1.  This  algorithm  consists 
of  transforming  between  object  and  Fourier  domains  and  applying 
appropriate  constraints  in  the  respective  domains.  He  use  the  symbol 
gi^(x)  to  represent  the  estimate  of  a  object  given  by  the  kth  iteration 
of  the  ER  algorithm.  The  prime  notation  In  9|^'(x)  Indicates  a  version 
of  the  kth  estimate  for  which  the  Fourler-domain  constraints  have  been 
enforced.  He  use  uppercase  symbols  to  denote  a  Fourler-domam 
representation  of  a  function.  In  practice  the  data  are  always  sampled 
and  therefore  we  use  the  discrete  Fourier  transform  (DFT) 


6(u)  • 


(6-1  ) 


and  Its  inverse 


9(«)  •  Giu).’'""”''"  6-: 

a 

in  the  algorithm.  Of  course  the  OFT  is  most  efficiently  computed  Mith  a 
fast  Fourier  transform  (FFT).  In  Eqs.  (6-1)  and  (6-2)  x  and  u  are 
two-dimensional  vectors  in  the  object  ano  Fourier  domains,  respectively, 
and  the  summation  notation  is  understood  to  represent  a  separate 
suMMtion  for  each  component  of  the  vector  running  from  0  to  N-1. 

In  order  to  '*nforce  a  given  constraint  we  define  a  least-squares 
constraint  operator.  The  function  of  the  operator  is  to  produce  an 
output  that  conforms  to  the  constraint  but  differs  from  the  input  as 
little  as  possible  in  a  least-squares  sense.  It  can  be  shown  that  when 
the  constraint  operators  have  this  property  the  mean  squared  error 
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between  the  latest  estimate  and  the  data  or  known  Information  will 
decrease  (or  stay  the  same)  at  each  Iteration  [6.3].  Thus  as  the 
algorithm  proceeds,  the  reconstruction  estimate  conforms  more  and  more 
closely  to  the  given  constraints.  This  Is  the  motivation  for  the  title 
"error  reduction." 

Although  many  types  of  constraints  have  been  used  with  the  ER 
algorithm,  our  problem  affords  a  modulus  constraint  In  the  Fourier 
domain  In  conjunction  with  a  support  constraint  In  the  object  domain. 
This  specific  realization  of  the  ER  algorithm  is  Illustrated  In  Figure 
6-2.  The  modulus  constraint  Is  performed  by  substituting  the  modulus  of 
the  latest  estimate  with  the  known  modulus  while  leaving  the  Fourier 
phase  untouched: 


G-(u)  -  (6-3) 
iG(u)| 

where  F(u)  Is  the  known  Fourier  modulus.  The  object  domain  constraint 
IS  equally  straightforward  and  is  enforced  by  setting  the  values  of  all 
pixels  that  fall  outside  of  the  support  equal  to  zero: 


9^' (x)  ,  X  £  S 

0  .  xe  S’ 


(6-4) 


where  S  stands  for  the  set  of  pixels  within  the  known  support  and  S'  is 
the  complement  of  S. 

In  order  to  monitor  the  progress  of  the  ER  algorithm  It  is  useful 
to  define  an  error  metric  for  each  of  the  constraints.  The  error  metric 
Is  essentially  a  mean  squared  error  between  estimates  before  and  after  a 
constraint  has  been  applied  and  Indicates  the  degree  of  agreement 
between  the  latest  estimate  and  the  known  constraint.  The  error  metric 
for  the  Fourier  modulus  constraint  Is  defined  as  follows: 
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FIGURE  f-2.  ERROR  REDUCTION  ALGORITHM  FOR  FOURIER  MODULUS  AND  OBJECT 
SUPPORT  CONSTRAINTS 
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(6-5) 


ep^  •  N'^  [;G(u)|  -  |F(u)|]^ 

U 

The  error  metric  for  the  object-support  constraint  Is  given  by 

*0^  -  S  I  g’(x)|^  (6-6) 

xeS'  '  ' 

As  the  algorithm  proceeds  both  of  these  error  metrics  Mill  decrease.  If 
they  simultaneously  achieve  values  close  to  or  equal  to  zero  then  the 
algorithm  has  achieved  a  restoration  that  has  good  agreement  with  both 
constraints. 

2 

Suppose  that  we  treat  the  error  metric  e^  as  an  objective  function 

to  be  used  in  a  gradient-search  algorithm.  Our  desire  is  to  minimize 

the  objective  function  by  varying  a  set  of  parameters  in  the  estimate. 

The  parameters  we  employ  are  the  individual  pixel  values  of  the 

estimate.  For  the  present  we  treat  only  real-valued  objects  which 
2 

require  N  independent  parameters  for  an  NxN  image  (complex  objects 

require  2N  parameters).  The  j  pixel  in  the  object  domain  is  located 

by  a  vector  x.  where  the  subscript  j  represents  any  convenient  ordering 
2  J  2 

of  the  N  pixels.  We  construct  an  N  -dimensional  Euclidian  vector 
space  for  which  each  coordinate  axis  corresponds  to  an  individual 
pa'^ameter.  Each  point  in  this  parameter  space  therefore  corresponds  to 
an  object  estimate  and  may  be  represented  by  the  parameter  vector  g(x). 
We  represent  the  j^^  parameter  and  its  associated  parameter  space  unit 
vector  by  9(Xj)  and  Vj,  respectively.  The  unit  vector  Vj  may  be 
interpreted  as  an  estimate  for  which  all  pixels  are  zero  except  for  the 
j^^  pixel  which  has  unit  strength.  This  vector  may  also  be  represented 
by  the  Kronecker  delta  6  .  The  objective  function,  ec  (g(x)),  1s  a 

A  X  9  A  J  r 

function  of  the  N  parameters,  and  may  be  visualized  as  a  surface  in  an 
2 

N  -t-l  -dimensional  space.  If  we  were  able  to  calculate  the  gradient  of 
this  surface  at  given  estimate  locations  then  well-known  gradient-search 
methods  could  be  employed.  The  gradient  is  formally  expressed 
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(6-7) 


Vep(g(x)) 


One  method  of  computing  the  gradient  Is  to  proceed  numerically  using  a 
finite  differences  approximation  to  the  partial  derivative: 

3ep^  ep^(g(x)  +  av.)  -  ep^(g(x)) 

- i — ^ 

where  a  Is  small  compared  with  significant  feature  sizes  in  the 

objective  surface.  This  brute-force  approach  Is  computationally 

2 

prohibitive  since  each  evaluation  of  ec  Involves  an  NxN  FFT  and  this 

2 

must  be  accomplished  for  each  of  the  N  parameters.  Fortunately  Flenup 
[6.3]  showed  that  the  exact  partial  derivative  may  be  calculated 
analytically  as  follows: 


We  reemphasize  that  the  prime  Indicates  that  the  Fourier  magnitude 
constraint  has  been  applied  to  the  estimate.  If  Eq.  (6-9)  is 
substituted  Into  Eq.  (6-7)  the  result  Implies  that  the  entire  gradient 
may  be  evaluated  with  a  forward  and  an  Inverse  FFT: 


Ve2(g(x))  .  21  2[g(Xj.)  -  g’(Xj)]vj 

*  2  [g(x)  -  g'(x)J 


(6-10) 


This  desirable  result  means  that  a  gradient  search  method  could 
realistically  be  employed  for  the  ep^(g(x))  objective  function. 
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Perhaps  the  simplest  gradient-search  algorltha  Is  the  Method  of 
steepest  descent  [6.10].  According  to  this  approach  the  latest  estimate 
May  be  Improved  upon  by  moving  In  parameter  space  In  a  direction 
opposite  that  of  the  gradient.  The  location  of  the  minimum  of  the 
objective  function  along  the  resulting  one-dimensional  cut  Is  then 
determined  giving  an  Improved  estimate.  This  pracedure  Is  repeated 
Iteratively  until  a  local  minimum  In  the  objective  function  Is  achieved. 

Some  optimization  problems  afford  additional  a  priori  Information 
about  disallowed  regions  In  parameter  space.  There  are  many  ways  of 
constraining  the  final  solution  to  the  allowed  region  of  parameter 
space.  One  obvious  way  of  Incorporating  this  Information  Is  to  proceed 
as  usual  with  the  steepest-descent  algorithm  until  an  estimate  Is 
produced  that  violates  the  a  priori  knowledge.  A  constraint  operator  Is 
then  employed  to  find  the  closest  allowed  estimate.  The 
steepest-descent  algorithm  Is  then  applied  to  the  latest  allowed 
estimate.  Unfortunately  this  constrained  steepest-descent  algorithm  can 
be  very  slow  since  the  direction  of  steepest  descent  Is  often  In 
competition  with  the  direction  enforced  by  the  constraint  operator. 

A  careful  analysis  of  the  ER  algorithm  reveals  that  It  Is,  In  fact, 
a  constrained  steepest-descent  algorithm  for  which  the  objective 
function  Is  ep  (g(x))  and  the  knowledge  of  object  support  defines  a 
disallowed  region  In  parameter  space.  The  kth  Iteration  of  the  ER 
algorithm  begins  with  an  estimate,  g|^(x),  and  replaces  Its  Fourier 
modulus  with  the  known  Fourier  modulus  to  get  g|('(x).  Notice  that  this 
Intermediate  result  Is  equivalent  to  moving  from  In  parameter 
space  In  a  9|('(x)  •  g|^(x)  direction;  that  Is,  1r  a  direction  opposite 
to  that  of  the  gradient.  In  fact  It  can  be  shown  that  the  objective 
function  Is  a  minimum  (zero)  at  g|^'(x).  Typically  9|('(x)  will  violate 
the  known  support  and  therefore  exists  In  a  disallowed  region  In 
parameter  space.  Applying  the  support  constraint  to  9|('(x)  produces  a 
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new  estimate,  i  that  now  resides  In  the  allowed  region,  thus 
completing  one  Iteration  of  the  constrained  steepest-descent  algorithm. 

While  we  have  thus  far  treated  the  Fourler-domaln  error  metric  as 
an  objective  function  we  could  Just  as  easily  have  selected  the  object- 
domain  error  metric,  e^  (g'(x)),  for  that  role.  The  gradient  for  this 
objective  function  Is  easily  obtained  because  the  calculation  of  the 
partial  derivative  with  respect  to  a  pixel  value  Is  more  direct: 


3gtXj) 


3 

sgtxj) 


fo  Xj  c  S 

[zgtxj)  xj  E  S’ 


(6-11) 


Recall  that  the  support  constraint  operator  sets  to  zero  all  pixels  in 
S'  and  leaves  those  In  S  untouched.  Clearly,  this  operation  moves  the 
latest  estimate  g|^'(x)  «  direction  opposite  that  of  Ve^  (g'(x)).  In 

addition  this  objective  function  Is  quadratic  along  this  one-dimensional 
cut  with  a  minimum  value  (zero)  at  gjj^^(x).  The  Fourier  modulus 
constraint  may  now  be  Interpreted  as  the  operator  that  takes 
out  of  a  new  disallowed  region  In  parameter  space.  Thus  the  ER 
algorithm  qualifies  as  a  constrained  steepest-descent  algorithm  from 
this  new  perspective  as  well. 

6.3  THE  SUMMED  OBJECTIVE  FUNCTION 

Historically  the  error  metric  e^  has  been  used  to  evaluate  an 
estimate  for  which  the  Fourier  constraints  have  been  satisfied. 
Consequently,  this  error  metric  is  a  function  of  the  pixel  values  In  a 
primed  estimate,  as  defined  In  Eq.  (6-6).  A  simple  generalization  of 
this  definition  yields  a  new  error  metric  that  can  be  applied  to  any 
estimate  g(x): 
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(6-12) 


ef(g(x))-  S  Cg(x)]^ 
0  xeS' 


2 

It  1s  easy  to  show  that  the  partial  derivative  of  with  respect  to 

pixel  values  In  the  estimate  has  the  same  form  as  given  In  Eq.  (6-11). 

Clearly  this  generalized  objective  function  and  Its  gradient  still 

pertain  to  functions  for  which  the  Fourier  constraints  have  been 

satisfied.  Notice,  however,  that  e.  (g(x))  now  has  the  same  underlying 
2  ° 

parameters  as  ep  (g(x)).  This  observation  affords  still  a  third 
Interpretation  of  the  ER  algorithm  that  yields  new  Insight.  The  ER 
algorithm  may  be  viewed  as  alternately  performing  steepest-descent 
operations  on  two  obJect1v<i  functions*  ep  (g(x))  and  (g(x)),  that 
coexist  In  the  same  parameter  space.  In  practice  It  Is  often  observed 
that  the  ER  algorithm  converges  rapidly  for  Iterations  early  In  the 
sequence  but  that  convergence  becomes  painfully  slow  as  the  Iteration 
number  Increases.  This  Is  because  the  work  performed  In  minimizing  the 
e^g(x})  objective  function  Is  largely  nullified  when  minimizing  the 
e^(g(x))  objective  function,  and  vice  versa.  Figure  6-3a  Illustrates 
this  point  pictorlally.  This  viewpoint  suggests  the  definition  of  a  new 
objective  function  that  Is  the  sum  of  the  opposing  objective  functions: 

e5^(g(x))  »  ep^(g(x))  ♦  e^^(g(x))  .  (6-13) 

The  gradient  of  this  new  objective  function  Is  simply  the  sum  of  the 
gradients  already  derived: 

^5^(g(x))  ■  7ep^(g(x))  +  7eQ^(g(x))  .  (6-14) 

The  calculation  of  this  new  gradient  Involves  a  forward  and  Inverse  FFT 
and  a  small  amount  of  computational  overhead.  Figure  6-3b  suggests  how 
moving  In  a  direction  opposite  that  of  the  gradient  of  the  sunwed 
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FIGURE  6-3.  OBJECTIVE  FUNCTION  SURFACES  FOR  TWO  PARAMETER  OBJECTS 
a.  Surfaces  used  in  error  reduction,  b.  Sunned  objective  function 
surface. 
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objective  function  may  circumvent  stagnation  due  to  opposing 

constraints.  Notice  that  If  we  choose  to  remain  with  steepest  descent 

2 

using  e^  ,  the  stepsize  still  has  to  be  determined.  This  can  be 
accomplished  by  one  of  a  variety  of  line  search  methods  that  utilize 
additional  samples  of  the  objective  function.  Each  additional 
objective-function  evaluation  requires  a  single  forward  FFT. 
Furthermore,  because  the  gradient  of  the  summed  objective  function  Is  so 
easily  computed  more  sophisticated  gradient-search  methods  such  as  the 
method  of  conjugate  gradients  or  a  memoryless  quasi-Newton  method  [6.10] 
may  profitably  be  employed.  Finally,  a  simple  generalization  of  these 
Ideas  to  Include  complex  objects  Is  found  In  Appendix  E. 

6.4  THE  eQ2(g(x))  OBJECTIVE  FUNCTION 

We  now  briefly  review  the  basic  characteristics  of  the  so-called 
Input-output  phase-retrieval  algorithms.  These  observations  will 
suggest  the  defining  of  a  new  objective  function  that  will  serve  as  an 
alternative  to  the  summed  objective  function. 

It  Is  convenient  to  partition  an  Iteration  of  the  ER  algorithm  Into 
two  steps.  The  first  step  enforces  the  Fourler-domain  constraints  while 
the  second  step  enforces  the  object-domain  constraints.  For  the  moment 
we  focus  on  the  first  step.  This  step  Involves  a  Fourier  transformation 
of  the  latest  estimate,  a  substitution  of  the  Fourier  modulus  by  the 
known  values,  and  an  Inverse  Fourier  transformation.  Together,  these 
operations  constitute  the  enforcement  of  Fourier  knowledge  and  may  be 
viewed  as  a  single  nonlinear  operation.  This  Is  depicted  schematically 
In  Figure  6-4.  It  Is  Important  to  recognize  that  any  output  of  this 
operation  will  satisfy  the  Fourler-domain  constraints  and  consequently 
ep  will  be  zero.  Should  the  output  also  satisfy  the  object-domain 
constraints  then  a  solution  has  been  found.  This  suggests  that  clever 
adjustments  to  the  Input  function  might  produce  an  output  that  more 

closely  satisfies  the  object-domain  constraints.  The  degree  of 

2 

consistency  with  the  support  constraint  can  be  monitored  by  the  e^ 
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error  metric  defined  In  Eq.  (6-6).  A  variety  of  feedback  strategies 

borrowed  from  nonlinear-systems  control  theory  can  be  employed  to  modify 

2 

the  latest  Input  In  order  to  drive  the  e^  error  metric  toward  zero. 
The  use  of  each  feedback  rule  defines  an  Individual  algorithm  and  the 
collection  of  feedback  rules  defines  the  class  of  Input-output 
phase-retrieval  algorithms.  All  feedback  rules  that  have  been  employed 
to  date  are  point  operations  meaning  that  an  Input  pixel-value 
adjustment  Is  based  solely  upon  the  desired  change  In  the  corresponding 
output  pixel  value. 

We  recognize  Immediately  that  If  a  solution  were  to  serve  as  an 
Input  function  then  It  will  pass  through  the  nonlinear  modulus  operator 
unchanged.  Notice  however  that  other  Inputs  can  also  output  a  solution. 
In  fact  any  Input  function  with  the  proper  Fourier  phase  will  produce  a 
solution.  Thus  a  solution  will  result  from  any  of  an  uncountable 
Infinity  of  Input  functions,  many  of  which  differ  dramatically  from  the 
solution.  The  ER  algorithm  may  be  viewed  as  a  particular  Input-output 
algorithm  for  which  the  feedback  rule  drives  the  Input  (as  well  as  the 
output)  toward  a  solution.  By  contrast  most  Input-output  algorithms 
have  a  more  flexible  feedback  rule  since  they  may  converge  upon  any  of 
the  many  Input  functions  that  yield  a  true  solution  upon  output. 

2 

He  reiterate  that  any  output  for  which  e^  Is  zero  will  be  a 
solution.  Therefore,  the  task  of  simultaneously  minimizing  the  Fourier 
and  object-domain  error  metrics  has  been  converted  Into  minimizing  a 
single  error  metric.  Unlike  the  summed  objective  function,  however, 
this  blending  of  the  two  error  metrics  Into  one  is  accomplished  without 
resorting  to  ad-hoc  methods  such  as  summing. 

We  use  the  term  objective  function  to  refer  to  an  error  metric  in 
conjunction  with  a  set  of  underlying  parameters.  A  logical  candidate 
for  an  alternative  objective  function  suggested  by  Input-output 
algorithms  Is  the  e^  error  metric  as  a  function  of  Input  pixel  values. 
This  new  objective  function  should  not  be  confused  with  the  e^^(g'(x)) 
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2 

Objective  function  used  In  the  ER  algorithm  which  treats  the  N 
object-estimate  (output)  pixel  values  as  parameters.  By  contrast  the 
new  objective  function,  e^  (g(x)),  utilizes  the  Input-function  pixel 
values  as  parameters  associated  with  the  object  estimate  given  upon 
output.  Having  made  this  subtle  but  critical  distinction  we  may  now 
write  an  expression  for  the  gradient  of  the  (g(x))  objective 
function: 


3e 


(6-15) 


As  before  a  numr/'ical  computation  of  the  gradient  Is  overwhelming.  It 
Is  natural  to  ask  If  an  analytic  expression  for  the  gradient  can  be 
derived.  While  the  details  of  this  calculation  are  outlined  In  Appendix 
F,  we  give  the  suprisingly  simple  result  here: 


J!L 

39(Xj) 


E 


|F(u)|Gj(u) 


6'(u)G*  (u)' 
G*(u) 


g12uu»Xj/N 


(6-16) 


where  *  denotes  complex  conjugate  end  6^(u)  Is  the  Fourier  transform  of 
an  error  Image  gj(x),  where 


g^(x)  ■  S'(x)g'(x) 


(6-17) 


and 


S'(x)  - 


1  ,  x  •  S' 
0  ,  X  E  S 


(6-18) 
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Three  FFT  operations  are  required  to  compute  G  from  g.  A  very 

6 

Important  feature  of  the  analytic  partial  derivative  quoted  In  Eq. 
(6-16)  Is  that  It  has  the  form  of  a  DFT.  The  Implication  Is  that  given 
the  expression  within  the  brackets  all  partial  derivatives  needed  to 
compute  the  gradient  are  provided  by  a  single  DFT.  Thus  the  total 
computational  cost  of  finding  7*^  (g(x))  for  a  given  Input  function  is 
four  FFTs  plus  minor  overhead.  With  these  manageable  computational 
requirements  the  e^  (g(x))  objective  function  may  be  minimized  via 
various  gradient-search  algorithms.  Some  care  must  be  taken  in  the 
evaluation  of  Eq.  (6-16)  to  avoid  division  by  zero.  This  problem  can  be 
circumvented  by  adding  a  small  constant  to  the  Fourier  magnitude  of  the 
Input  function  at  those  spatial  frequencies  for  which  |6(u)|  Is 
Identically  zero. 

Notice  that,  like  Input-output  algorithms,  there  are  many  Input 
functions  to  which  a  gradient-search  algorithm  can  converge  for  this 
objective  function.  This  means  that  the  objective  function  contains 
many  global  minima,  each  equally  acceptable  for  producing  a  solution  as 
an  output.  It  Is  conceivable  that  this  multiplicity  of  Input  solutions 
could  yield  faster  convergence  rates  than  an  objective  function  having 
only  a  single  global  minimum  (e.g.  the  summed  objective  function). 

It  Is  useful  to  recognize  that  any  gradient-search  algorithm  used 
In  conjunction  with  the  (g(x))  objective  function  may  also  be 
Interpreted  as  a  particular  feedback  rule  In  an  Input-output  algorithm. 
Unlike  other  feedback  rules,  however,  this  rule  Is  not  a  point 
operation.  In  other  words,  the  gradient-search  feedback  rule  is  more 
flexible  than  other  existing  rules  since  many  Input  pixels  may  be 
adjusted  In  order  to  effect  a  desired  change  In  a  single  output  pixel. 

Unfortunately,  there  Is  no  guarantee  a  priori  that  the 
objective  function  has  a  surface  contour  that  lends  Itself  to 
minimization  via  gradient  search.  For  example  the  e^^(g(x))  surface 
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FIGURE  6-5.  PRELIMINARY  IMAGES  DERIVED  FROM  MINIMIZING  THE  e 
OBJECTIVE  FUNCTION 


My  contain  Mny  local  MlniM  In  which  gradlant-saarch  algorlthiis  could 
bocoM  tntrapptd.  Answtrs  to  such  questions  art  often  the  byproduct  of 
extensive  experlMntatlon. 

2 

Soiie  prellalnary  experlMnts  were  perforMd  In  which  the  e^  (g(x)) 
objective  function  was  used  In  conjunction  with  the  Mthod  of  steepest 
descent.  A  nuaber  of  observations  can  be  eade  about  the  results 
displayed  In  Figure  6-5.  Notice  the  doalnant  stripes  In  the  gradient 
iMge  for  the  first  Iteration.  By  gradient  iMge  we  Man  the  iMge  for 
which  each  pixel  value  Is  assigned  the  value  of  the  associated  coeponent 
of  the  gradient.  This  Is  the  Inage  that  Is  scaled  and  added  to  the 
latest  Input  laage  to  acquire  the  succeeding  input  Inage  In  a 
steepest-descent  scheM.  These  stripes  are  Intriguing;  but  their  origin 
Is  unknown  at  present.  The  aagnltude  of  the  gradient  was  observed  to 
decrease  with  Iteration  nunber.  As  a  result,  the  stripes  fron  the  first 
gradient  iMge  still  persist  In  the  100th  Input  iMge.  Notice,  however, 
that  the  stripes  do  not  appear  In  an  output  iMge,  which  Is  consistent 
with  the  notUn  that  the  Input  iMge  need  not  resenble  the  output  Image. 
It  Is  encouraging  that  after  100  Iterations  the  output  iMge  bears  a 
rough  resemblance  to  the  true  object.  Nore  experlMntatlon  with  this 
objective  function  Is  needed  before  a  judgement  can  be  made  about  Us 
usefulness.  For  example,  more  sophisticated  gradient-search  Mthods 
would  have  a  better  chance  of  converging  to  a  solution.  Should  the 
Cq  (g(x))  objective  function  In  conjunction  with  the  best 
gradient-search  methods  prove  not  to  be  competitive  with  current 
input-output  algorithms.  It  may  yet  be  useful  for  breaking  out  of 
stagnation  episodes. 

He  conclude  this  section  by  noting  that  while  we  have  restricted 
objects  to  be  real-valued  for  simplicity,  the  case  admitting  complex 
objectr.  Is  of  great  Interest  when  objects  are  Illuminated  coherently. 
The  definition  and  derivation  of  the  gradient  of  the  objective 
function  for  complex  objects  Is  presented  In  Appendix  6. 
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FIGURE  6-5.  PRELIMINARY  IMAGES  DERIVED  FROM  MINIMIZING  THE 
OBJECTIVE  FUNCTION 
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6.5  FOURIER  PHASE  PARAMETERS 

The  choice  of  underlying  parameters  for  an  objective  function  can 
have  a  tremendous  Impact  upon  the  behavior  of  gradient-search 
algorithms.  To  this  point  we  have  selected  the  Input  pixel  values  (or 

real  and  Imaginary  parts  of  the  Input  pixels)  as  our  N  (or  2N  ) 

2 

parameters  underlying  the  e^  (g(x))  objective  function.  This  choice  has 
merit  since  It  affords  an  analytic  expression  for  the  gradient  requiring 
only  four  FFTs.  An  alternative  and  very  different  set  of  parameters 
worth  consideration  Is  the  set  of  Fourier  phase  values  In  a  Fourier 
estimate  of  a  solution.  Because  the  Fourier  modulus  is  known,  a  Fourier 
estimate  Is  determined  by  an  estimate  of  the  Fourier  phase,  ')>(u): 

G'(u)  •  |F(u)|«'*<“>.  (6-19) 

An  Inverse  FFT  gives  the  corresponding  object-domain  estimate, 

g'(x)  - 

u 

This  estimate  may  also  be  Interpreted  as  the  output  from  an  input-output 
algorithm  since  It  has  the  proper  Fourier  modulus.  Consequently,  the 
object-domain  error  metric  can  be  computed: 

*0  *  |g’(x)!^  (6-21) 

xeS' 

2 

The  e^*^  error  metric  1s  therefore  implicitly  a  function  of  the  Fourier 

phase  values  and  e^^((^(u))  serves  as  the  third  new  objective  function 

Introduced  In  this  chapter.  We  mention  parenthetically  that  throughout 

this  section  we  allow  for  complex-valued  objects  since  there  Is  no 

simplification  of  derivations  by  resorting  to  real-valued  objects. 

Notice  that  the  designation  of  the  Fourier  phase  values  as  the 

2 

underlying  parameters  has  fixed  the  number  of  parameters  at  N  .  This 
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Is  exactly  half  the  number  of  parameters  that  occur  when  using  the  real 
and  Imaginary  parts  of  the  Input  pixel  values  as  parameters.  It  remains 
to  be  seen,  though.  If  an  analytic  expression  for  the  gradient  of  the 
object-domain  error  metric  with  respect  to  the  Fourier  phase  parameters 
can  be  derived. 


The  gradient  Is  defined  as 


Ve2((j)(u)) 


y  f^O _ 

^  ''j 


(6-22) 


where  is  the  unit  vector  In  parameter  space  associated  with  the 
phase  parameter  at  location  Uj  In  the  Fourler-domain  estimate.  As 
usual,  the  heart  of  the  gradient  Is  the  partial  derivative 


3<)»(Uj) 


3 

3(Ji(Uj) 


xeS* 


(6-23) 


(6-24) 


where  C.C.  stands  for  complex  conjugate.  The  partial  derivative  In  Eq. 
(6-24)  may  be  simplified: 


3g'*(x) 

94>(Uj) 


U 


3  g-4(u) 

3<l»(Uj) 


,-2 


E  |p(u)i 


g-12Tru*x/N^_^  j^-l^(u)^ 


u,u. 


(6-25) 

(6-26) 
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Applying  the  sifting  property  of  the  Kronecker  delta, 
(6-26)  leaves  only  one  tern  froei  the  sumaatlon: 


In  Eq. 


(6-27) 


Substituting  back  Into  Eq.  (6-24): 


(6-28) 


j^g'(x)(-i)fr^|F(Uj)|e'^^’^“j’*/^e’’^(‘^j)  +  c.C.j 

=  N*2  |F(Uj.)|jj^(-i)e"^^^‘'j^  Z  S'(x)g'(x)e"’^''‘'j*’‘^^^  +  C.C.j 


(6.29) 


The  summation  In  Eq.  (6-29)  Is  the  Fourier  error  Image,  6g(u),  defined 
In  the  previous  section  by  Eq.  (6-17).  Therefore  we  have 


3  ^ 

-  -  N-2  |F(uJ|r(-1)G.(uJe"^^^“j^  +  C.C.l  (6-30) 

30(Uj)  J  L  «  J  J 

>  2N‘2  |F(Uj)|Im  (6-31) 
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where  Im|*| stands  for  the  Imaginary  part. 

Again  we  have  been  able  to  find  an  expression  for  the  gradient  with 
a  remarkably  compact  form.  Equation  (6-31)  Implies  that  the  component 
of  the  gradient  associated  with  the  spatial  frequency  Uj  Is 
proportional  to  the  modulus  at  that  spatial  frequency  and  Is  dependent 
upon  the  Fourler-domain  error  Image  and  the  latest  Fourier  phase  In  a 
less  direct  way.  An  examination  of  Eq.  (6-31)  reveals  that  the  entire 
gradient  can  be  computed  with  2  FFTs  plus  minor  overhead.  The  actual 
evaluation  of  the  objective  function  for  a  particular  Fourler-phase 
estimate  requires  only  one  FFT.  Thus  employing  Fourler-phase  values  as 
optimization  parameters  Is  certainly  competitive  with  the  use  of 
Input-pixel  values  from  the  standpoint  of  operations  required  to  compute 
the  gradient.  How  these  two  gradient-search  formulations  compare  with 
respect  to  convergence  properties  can  only  be  determined  by 
experimentation.  We  might  expect  the  Fourler-phase  formulation  to 
perform  differently  since  the  Fourler-phase  parameters  are  so  different 
In  character  from  and  nonllnearly  related  to  the  Input  pixel-value 
parameters.  Use  of  the  Fourier  phase  for  parameters  has  the  added 
appeal  that  these  are  In  fact  the  unknowns  In  the  phase-retrieval 
problem.  As  a  result  the  Fourier  phase  formulation  Is  somewhat  mure 
direct  and  may  lend  Itself  to  analysis  when  noise  Is  present. 

6.6  CONCLUSIONS  AND  FUTURE  WORK 

We  have  shown  that  the  ER  algorithm  may  be  Interpreted  as  a 
constrained  steepest-descent  algorithm  for  which  the  objective  function 
consists  of  the  Fourler-domain  error  metric  as  a  function  of  pixel 
values  In  the  latest  estimate.  In  addition  we  have  proposed  three  new 
objective  functions  for  performing  phase  retrieval  using  gradient-search 
methods.  These  Include:  1)  use  of  the  summation  objective  function  with 
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pixel  values  of  the  latest  estimate  as  optimization  parameters,  2)  use 
of  the  object-domain  error  metric  with  Input  pixel  values  as  parameters, 
and  3)  use  of  the  object-domain  error  metric  with  Fourler-phase  values 
as  parameters.  Analytic  expressions  for  the  gradients  for  each  of  these 
approaches  have  been  derived.  The  simplicity  of  these  expressions 
Implies  that  gradient-search  methods  have  the  hope  of  being 
computationally  tractable  and  even  competitive  with  existing 
Input-output  algorithms.  The  total  number  of  FFTs  required  to  evaluate 
the  objective  function  and  compute  the  gradient  for  each  of  these 
approaches  Is  shown  In  Table  6.1. 

TABLE  6.1.  NUMBER  OF  FFTs  REQUIRED  FOR  GRAD  I  ENT- SEARCH  APPROACHES 

Objective  Function  IFFTs  to  evaluate  IFFTs  to  evaluate 

objective  function  gradient 

e|(g(«))  I  2 

«,^8(«))  2  5 

I  2 

Of  course  extensive  experimentation  needs  to  occur  to  see  if  the 
surface  contour  of  each  proposed  objective  function  Is  well  suited  for 
gradient-search  methods.  Surface  contour  depends  upon  such  things  as 
the  Intrinsic  definition  of  the  objective  function,  the  particular  true 
object,  and  the  amount  of  noise  In  the  data.  The  suitability  of  a 
particular  gradient-search  algorithm  to  a  given  surface  contour 
manifests  Itself  In  the  convergence  rates  of  the  algorithm.  For  example 
a  memory  less  modified  Newton  method  [6.10]  may  converge  well  with  the 
same  objective  function  for  which  a  steepest-descent  algorithm 
stagnates.  Software  Is  currently  being  developed  to  test  for  the 
convergence  rates  as  a  function  of  type  of  objective  function,  choice  of 
true  object,  amount  of  noise  and  gradient-search  method  employed.  In 
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addition,  th«st  gradltnt-starch  approachts  nttd  to  ba  tastad  in  a  rola 
that  coHplaaants  currant  Input-output  algorlthas.  Gradlant-saarch 
approachas  could  Mka  a  significant  contribution  to  tha  flald  of  phasa 
I'atrlaval,  should  thay  consistantly  provida  a  aoda  of  ascapa  fron  any  of 
tha  various  typasof  stagnation  that  hava  baan  known  to  appaar  with 
Input-output  algorlthas. 
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7 

MODELING  APPROACH  TO  PHASE  RETRIEVAL 

Tht  aodtllng  approach  Is  a  ntw  Mthod  for  atttnptlng  to  solvt  the 
phast  rttrleval  problan.  In  this  sactlon  wo  dtscrlba  tht  nodal Ing 
approach  In  gtntral  tarns,  and  than  discuss  a  particular  lapleaentatlon 
that  was  attaaptad. 

Lot  F(u,v)  ■  |F(u,v)|  axpCU(u,v)]  ba  tha  coaplax  Fourlar  transform 
of  a  particular  objact.  Supposa  that  althar  F(u,v)  ovar  tha  antlra 
maasuraaant  apartura  or  F(u«v)  ovar  soaa  saall  arta  can  ba  aodeled  by  a 
paraaatcrixad  function.  M: 

N(u.v,;a.b,...)  -  lN(u,v;a.b...)|  axp[1^(u,v;a.b,...)],  (7-1) 

wharf  a«b«...  art  unknown  paraaatars.  If  wa  art  givan  only  tha  Fourlar 
aodulus*  |F(u.v)|,  than  It  night  ba  possibla  to  astinata  tha  phasa. 
l'(u,v),  by  (1)  finding  tha  valuas  of  tha  paranatars  a.b,...  that  bast 
fit  tha  modulus  of  tha  nodal,  |N(u,v;a,b,. . . ) |  ,  to  |F(u,v)|  ,  and  (2) 
evaluating  (;)(u,v;a,b,...)  for  that  sat  of  valuas  of  tha  paranatars. 

Tha  nost  difficult  part  of  this  approach  Is  finding  a  nodal,  M, 
that  Is  suitable. 

In  a  first  attanpt  at  using  tha  nodal Ing  approach,  each  snail  area 
about  tha  local  naxina  of  tha  Fourlar  modulus  was  nodalad  using  a 
function  taken  from  tha  control  theory  literature.  Supposa  that 
contours  about  a  local  naxinun  of  tha  Fourlar  nodulus,  at  a  level  3  dB 
down  from  tha  local  naxinun,  have  an  elliptical  shape,  with  tha  major 
axis  of  length  w^^  at  an  angle  relative  to  tha  u*ax1s  and  minor  axis 
w^.  Let  tha  local  naxinun  ba  at  location  (U|^,V|^)  where  It  has  tha  value 

•  I  •  (7-2) 
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Also  define  the  dlstence  fron  a  given  point  (u.v)  to  the  peak  as 

7  2 

«  -  [(u-Uj,)2  ♦  (v-v^)^]  .  (7-3) 


and  let 

8  ■  tan"^[{v-Vjj)/(u-U|j)]  , 


(7-4) 


and 


Wc  ■  Wjj  cos(ejj-8)  ♦  s1n(e|^-e)  . 

Then  the  aodel  we  used  for  a  region  about  the  local  auxlnun  Is 

A.w  ^ 

M(e;Aj^,0j^,Wj^,w^,O)  ■  — - - 2 - 

+  i2ww^D 

which  has  squared  modulus 


-  W^)^  +  (2ww^D)^ 


and  phase 


♦  ■  -tan"^[2ww^D/(w^^  -  w^)]  . 


(7-5) 


(7-6) 


(7-7) 


(7-8) 


Note  that  the  parameters  w^,  w^  and  are  contained  within  w^. 
These  expressions  were  used  In  the  following  way: 

(1)  A  local  maximum  of  the  squared  Fourier  modulus  was 
found. 

(2)  A  curve  fit  of  Eq.(7-7)  to  the  squared  Fourier  modulus 
was  performed  to  estimate  the  unknown  parameters. 
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(3)  Tht  phasf  In  that  region  Mas  conputed  by  Eq.(7-8)  using 
the  paranetcr  estlaatas. 

(4)  Eq.(7-7)  Mas  evaluated  using  the  paraMter  estleates 
and  subtracted  froa  the  squared  Fourier  noduluSf 
leaving  the  residual  Fourier  Modulus. 

(5)  Repeat  steps  (1)  to  (4)  replacing  the  squared  Fourier 
Modulus  Mith  the  residual  Fourier  Modulus,  until  all 
the  Major  local  Maxlna  are  accounted  for. 

(6)  Fona  the  net  Fourier  phase  as  the  sum  of  all  the  phase 
functions  obtained  In  step  (3). 

(7)  Fora  an  laage  by  Inverse  Fourier  transforalng  the 
coaplex  function  foraed  froa  the  given  Fourier  Modulus 
and  the  net  Fourier  phase. 

2 

Note  that  for  large  m.  |N|  approaches  zero  and  ^  In  Eq.(7-8) 
approaches  zero,  so  the  Model  has  strong  local  effect  near  each  local 
maxiMUM  and  a  Meaker  effect  on  neighboring  points. 

When  the  procedure  Mas  perforaed  for  a  SAR  laage  of  the  type  used  In 
the  digital  experlnents  described  In  Section  5,  the  reconstructed  Image 
bore  no  resemblance  to  the  original  object.  The  reason  for  failure  Is 
not  totally  understood,  but  mc  speculate  that  the  model,  Eq.(7-6),  Is 
not  appropriate  to  the  Fourier  transforms  of  SAR  Images. 

If  further  Mork  along  these  lines  Mere  to  be  pursued.  It  eould  be 
Important  to  first  develop  more  appropriate  models  for  SAR  signal 
histories. 


no 
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LABORATORY  EXPERIMENTS 


The  objective  of  Task  3  of  the  program  Is  to  perform  laboratory 
experiments  which  demonstrate  reduced  tolerance  Imaging.  These 
experiments  will  validate  the  theoretical  developments  concerning 
constraints*  measurements,  phase  retrieval  and  Image  reconstruction 
algorithms,  and  uniqueness  and  sensitivity  Issues  under  more  realistic 
conditions  than  Is  possible  In  the  computer  simulations  performed  under 
Tasks  1  and  2.  The  use  of  real  objects.  Illumination  sources,  optics, 
and  detectors  will  place  greater  demands  on  the  reconstruction 
algorithms.  The  quality  of  the  reconstructed  Images  from  experimented 
data  will  be  compared  both  to  Images  from  computer  simulations  of  the 
experiment  and  to  "ground  truth"  Images  collected  In  the  laboratory  with 
a  conventional  sensor  having  an  equivalent  aperture.  Experimental 
parameters  (e.g.,  measurement  s1gna1-to>no1 se  ratio,  type  of  shape 
constraint,  sharpness  of  shape  constraint)  will  be  varied  for 
comparlslon  to  theoretical  and  computer  simulation  results. 

Two  experiments  simulating  different  types  of  systems  will  be 
performed:  an  active  coherent  experiment  In  the  visible  and  a  passive 
Incoherent  experiment  In  the  visible  or  Infrared.  The  results  of  the 
active  coherent  experiment  will  be  useful  for  predictions  concerning  SAR 
and  active  laser  Imaging  sytems,  whereas  the  passive  Incoherent 
experiment  will  be  pertinent  to  conventional  passive  and  passive 
Interferometric  Imaging  systems.  In  the  active,  coherent  experiment  a 
target  will  be  Illuminated  with  a  laser  (with  various  Illumination 
shapes)  and  Intensity  data  will  be  collected  In  the  far-fleld  of  the 
target.  Reconstruction  algorithms  will  then  be  used  to  determine  the 
phase  In  the  far-fleld  of  the  target  and  therefore  an  Image  of  the 
target.  Thus,  this  experiment  simulates  an  active,  coherent 
reduced- tolerance  sensor  which  measures  Intensities  only.  The  optical 
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and  electronic  equipment  requirements  of  this  experiment  have  been 
determined  and  most  of  the  additional  equipment  has  been  purchased, 
delivered,  and  tested.  An  Initial  setup  has  been  made  of  part  of  the 
equipment  In  the  laboratory.  Computer  software  Is  under  development  to 
collect  and  process  the  measurement  data.  The  two  experiments  will  be 
performed  serially,  so  further  planning  of  the  passive  experiment  Is 
purposely  proceeding  slowly  until  the  active  experiment  Is  approximately 
one-half  complete.  Further  discussion  of  both  experiments  Is  given  In 
Sections  8.1  and  8.2. 

8.1  ACTIVE  EXPERIMENT 

The  objective  of  the  active  experiment  Is  to  demonstrate  Imaging  of 
a  coherently  Illuminated  target  from  Intensity-only  measurements  made  In 
the  far-fleld.  This  simulates  a  sensor  having  greatly  reduced  tolerance 
to  the  position  and  quality  of  Its  receiving  aperture  compared  to  a 
conventional  Imaging  sensor.  The  wide  range  of  parameters  which  must  be 
considered  In  planning  this  experiment  and  which  are  available  to  be 
varied  to  test  theoretical  developments  and  computer  simulation  results 
are  discussed  In  Section  8.1.1.  Many,  but  due  to  finite  resources,  not 
all,  of  the  parameters  discussed  will  be  exercised  In  the  actual 
experiments.  The  experiment  design  Including  both  optics  and  electronics 
Is  discussed  In  Section  8.1.?.. 

In  order  to  greatly  Increase  the  range  of  parameters  which  could  be 
Investigated  with  a  fixed  amount  of  manpower.  It  was  decided  at  the 
beginning  of  the  program  that  an  array  processor  would  be  purchased 
which  was  sufficiently  powerful  to  perform  the  Iterative  Fourier 
transform  algorithms  needed  for  Image  reconstruction  and  which  was 
compatible  with  one  of  ERIM's  laboratory  computers.  (The  use  of 
existing  VAX  facilities  would  have  required  time-consuming  transfer  of 
large  amounts  of  data  and  vastly  Increased  computational  costs.)  A 


112 


suitable  array  processor  was  ordered,  but  delays  In  Its  delivery  have 
forced  the  experimental  work  to  fall  behind  Its  Initial  timetable.  A 
less  powerful,  but  similar,  array  processor  has  been  obtained  on  loan 
from  the  manufacturer  (Mercury  Computer)  until  the  original  order  can  be 
shipped.  Current  work  on  the  laboratory  experiments  Is  concerned  with 
writing  software  to  collect  and  process  the  measurement  data,  further 
setup  of  optical  components,  and  final  equipment  purchases. 

8.1.1.  ACTIVE  EXPERIMENT  PARAMETERS 

A  wide  range  of  parameters  Is  available  to  be  varied  In  the  active 
experiment  to  test  theoretical  developments  and  computer  simulation 
results.  These  parameters  must  also  be  carefully  controlled  to  ensure 
meaningful  results.  The  most  Important  of  these  parameters  are 
discussed  In  this  section. 

The  pattern  of  Illumination  on  the  target  can  be  described  by  Its 
spatial  and  temporal  coherence,  shape,  sharpness  of  edges,  phase 
distribution,  angle  of  Incidence  on  the  target,  and  polarization.  All 
of  these  parameters  may  affect  the  quality  of  the  reconstructed  Image. 
Equally  Important,  they  may  take  on  different  values  depending  on  the 
application  being  simulated  In  the  laboratory.  In  an  application  where 
the  target  Is  actively  Illuminated  by  a  laser,  the  spatial  (transverse) 
and  temporal  (longitudinal)  coherence  lengths  may  be  less  than  the 
target  size.  The  laboratory  system  can  allow  Illumination  with  variable 
coherence  lengths  either  by  manipulating  the  spatial  coherence  of  a  gas 
laser  or  by  using  a  broader-band  dye  laser.  It  Is  known  from  ERIM 
Investigations  that  the  Illumination  pattern  shape  and  sharpness  of  the 
edges  affects  reconstruction  algorithm  convergence.  In  applications, 
the  range  of  Illumination  shapes  and  sharpness  of  ed^es  may  be  limited 
by  practical  considerations  on  the  transmitting  aperture.  The 
laboratory  system  must  accept  a  variety  of  masks  and  Image  them  onto  the 
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target  through  a  controllable  finite  aperture  to  control  Illumination 
shape  and  sharpness.  Since  the  Illumination  phase  may  not  be  constant 
over  the  target  In  practical  applications,  the  masks  will  need  to  be 
holographic  If  experimental  control  of  the  Illumination  pattern  phase 
distribution  Is  desired.  Applications  may  be  either  monostatic  or 
bistatic,  so  the  laboratory  system  should  allow  for  either.  The  target 
will.  In  most  cases,  partially  depolarize  the  Illumination.  The 
experimental  system  should  be  capable  of  making  measurements  of  the  two 
orthogonal  components  of  the  light  at  the  detector. 

The  target  parameters  Include  reflectivity  contrast  and  structure, 
surface  roughness,  surface  topography  (3-0  nature  of  target),  motion 
during  the  measurement  process,  and  noncoherent  background  Illumination. 
Practical  targets  will  vary  In  their  roughness,  although  nearly  all  will 
be  rough  at  visible  and  Infrared  wavelengths.  The  experimental  system 
should  primarily  jse  rough  targets  to  create  real  speckle  effects. 
However,  It  may  be  useful  to  use  smooth  targets  (film  transparencies  in 
a  liquid  gate)  In  setting  up  the  experiment  to  test  and  debug  the 
optical  and  electronic  components  and  software.  Real  targets  are 
three-dimensional,  but  to  varying  degrees.  A  variety  of  3-D  objects 
should  be  available  for  the  experiment.  A  practical  reduced  tolerance 
sensor  may  need  to  cope  with  target  motion  during  the  detector 
Integration  time.  This  effect  can  be  simulated  by  mounting  the  target 
on  an  motor-driven  translation  or  rotation  stage.  Any  real  target  will 
also  be  noncoherently  Illuminated  from  various  thermal  sources.  While 
this  Illumination  will  only  add  a  uniform  bias  to  the  far-fleld 
measurements.  It  should  be  Included  In  the  experimental  setup. 

The  propagation  path  between  the  target  and  the  sensor  can  be  of 
such  a  length  as  to  be  either  near  or  far  field  and  can  Include 
atmospheric  turbulence,  scattering  (aerosols,  fog,  smoke),  and 
absorption.  In  an  application,  any  or  all  of  these  effects  may  be 
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present.  The  optical  setup  of  the  experiment  can  allow  for  Insertion  or 
removal  of  lenses  to  give  either  near  or  far  field  conditions  at  the 
detector.  (It  must  be  noted  that  the  size  of  the  speckles  In  the 
measurement  plane  will  depend  on  the  sensor  distance.)  Turbulence  with 
the  proper  statistical  properties  Is  very  difficult  to  simulate  In  an 
Indoor  laboratory.  The  best  approach  (and  one  which  allows  reproducible 
turbulence)  Is  to  use  movable  phase  plates.  These  are  glass  plates  with 
controlled  thickness  variations.  Scattering  and  absorption  are  easier 
to  simulate  (e.g..  with  fog  chambers,  optical  narrow  band  filters). 
Since  their  main  effects  are  to  reduce  signal  levels  and  Increase 
detector  bias  light  levels,  they  can  also  be  studied  as  described  In  the 
next  paragraph  on  detector  parameters. 

The  most  Important  detector  parameters  are  signal  level,  type  of 
noise,  noise  level,  background  Illumination  level,  spatial  and  temporal 
sampling  rates,  polarization  detected,  nonllnearltles  In  response,  and 
nonuniformities  In  response,  bias  and  noise.  The  values  of  these 
parameters  are  crucial  to  the  viability  of  any  real  application.  To 
give  useful  results,  the  experimental  setup  must  be  able  to  vary  the 
signal  and  background  Illumination  levels,  simulate  various  sampling 
rates,  detect  orthogonal  polarizations,  and  create  nonuniform  background 
Illumination  levels.  The  type  of  noise,  noise  level,  nonuni forml ties  In 
detector  response  and  noise  (pattern  noise),  and  nonllnearltles  In 
response  will  be  primarily  determined  by  the  detector  chosen  for  the 
experiment.  An  Important  aspect  of  the  experiment  design  and  procedure 
should  be  to  measure,  calibrate,  and  correct  for  the  effects  of  these 
parameters  on  the  measured  data.  This  may  be  a  difficult  task  and  much 
practical  experience  will  be  gained  which  may  be  applied  to  other 
detectors  In  the  future. 
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Many  applications  Involving  active  llluMlnatlon  can  be  expected  to 
have  low  signal  levels.  To  adequately  simulate  these  applications,  an 
Image  Intensifler  should  be  used  before  the  detector.  Even  with  a 
thermal-nolse-l Imited  detector,  the  use  of  an  Image  Intensifler  may 
allow  operation  In  a  shot  (photon)  noise  limited  mode.  The  Image 
Intensifler  will,  of  course,  also  have  nonuniformities,  nonl  Inearl  ties, 
and  a  spatial  sampling  rate  which  must  be  measured  and  considered  In  the 
experimental  setup  and  data  analysis. 

The  effects  of  speckle  are  so  Important  to  an  active  coherent 
experiment  that  Its  parameters  deserve  to  be  discussed  separately.  The 
speckle  In  the  measurement  plane  will  have  Its  size  determined  chiefly 
by  the  size  of  the  Illuminated  region  on  the  target  and  by  the  optics 
(If  any)  placed  between  the  target  and  the  sensor.  Ideally,  the 
detector  must  sample  the  Intensity  speckle  pattern  at  the  Nyquist  rate 
(two  samples  per  speckle)  or  greater.  For  a  given  detector,  magnifying 
optics  may  be  necessary.  It  may  also  be  desired  to  Investigate  the 
effect  of  measurements  at  less  than  the  Nyquist  rate.  Some  speckle 
reduction  techniques  Include  averaging  the  Intensities  of  Images  from 
Independent  looks  (aspect  angles)  at  the  target.  The  experimental  setup 
should  be  capable  of  rotating  or  translating  the  target  In  discrete 
steps  to  generate  these  Independent  looks. 

8.1.2  ACTIVE  EXPERIMENT  DESIGN 

An  optical  setup  and  electronic  hardware  for  performing  the  active 
experiment  Is  shown  In  Fig.  8-1.  The  laser  source  Is  spatially 
filtered,  collimated,  and  used  to  Illuminate  a  masx  which  Is  Imaged  onto 
the  target  via  a  beamsplitter.  The  target  Is  In  the  front  focal  plane 
of  lens  and  the  far-fleld  distribution  of  the  light  from  the  target 
Is  obtained  In  the  back  focal  plane.  This  distribution  Is  Imaged  with 
magnification  by  lenses  end  (together  forming  an  afocal  telescope) 
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Figure  8-1.  (U)  Active  Experiment  Laboratory  Setup 


onto  the  detector.  Some  of  this  light  Is  redirected  by  a  beamsplitter 
or  a  removeable  mirror  to  lens  which  forms  an  Image  of  the  target  on 
a  second  detector.  Signal  processing  Is  used  to  digitize  the  detector 
signals,  perform  preprocessing  of  the  data,  and  to  perform  the  phase 
retrieval  and  Image  reconstruction  algorithms.  A  laboratory  computer  Is 
used  to  control  the  experiment  and  for  data  storage. 

The  experiment  design  outlined  above  permits  most  of  the  parameters 
discussed  In  Sec.  8.1.1  to  be  controlled.  For  example.  It  allows  the 
use  of  various  masks  and  targets,  target  rotation  and  translation, 
noncoherent  background  Illumination  of  target  and  detector,  variation  of 
detector  signal  levels,  matching  of  speckle  size  to  detector  resolution, 
detection  of  orthogonal  polarizations,  use  of  an  image  Intensifler,  and 
detection  of  the  target  Image  using  a  conventional  optical  sensor. 
Simple  rearrangements  of  the  optical  equipment  will  allow  holographic 
masks  to  be  used,  the  spatial  coherence  of  the  Illumination  to  be 
varied,  monostatic  or  bistatic  Illumination,  near-field  Intensities  to 
be  measured,  and  the  Inclusion  of  simulated  turbulence,  scattering,  and 
absorption. 

Target,  Optics,  and  Detector  Design 

For  the  active  (far-fleld)  experiment,  the  following  parameters 
need  to  be  chosen  or  determined: 

Target  diameter,  d 

Target  1-0  space-bandwidth  product,  N 
Lens  Lj  focal  length,  F^ 

Lens  L|  aperture  diameter,  0^ 

Magnification  by  lenses  L2  and  L^,  M 
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Far-fleld  detector  dlaaeter,  s 
Far-fleld  detector  element  spacing.  M 
Lens  L^  focal  length,  F^ 

Target  Image  detector  diameter,  w 
Target  Image  detector  element  spacing, 

Wavelength,  X. 

These  parameters  are  related  by  a  number  of  equations.  Which  parameters 
can  be  freely  chosen  and  which  are  thereby ' determi ned  depends  upon  one's 
point  of  view.  The  following  explanation  Is  based  on  the  fact  that,  for 
experiments  performed  In  this  program,  most  of  the  parameters  will  be 
determined  by  the  capabilities  of  available  detectors. 

The  allowable  target  space-bandwidth  product  (SBP)  is  obviously 
limited  by  the  number  of  detector  elements  In  the  far-fleld  detector. 
Because  of  the  squaring  operation  Inherent  In  Intensity  measurements, 
the  spatial  frequency  spectrum  of  the  detected  signal  Is  doubled 
(relative  to  amplitude  detection)  and  therefore  a  detector  with  2N 
elements  In  each  dimension  Is  required  to  collect  the  data  from  which  a 
(possibly  complex-valued)  target  Image  with  a  SBP  equal  to  N  (in  each 
dimension)  can  be  reconstructed.  For  current  solid  state  array 
detectors  operating  In  the  visible,  the  number  of  detector  elements  In 
either  dimension  varies  from  about  240  to  490.  Therefore,  assuming  that 
a  square  Image  with  N  equal  to  a  power  of  2  Is  desired,  N  will  be  no 
greater  than  128. 

If  the  target  field  diameter,  d.  Is  chosen,  then  the  resolution 
required  of  the  optical  system  to  the  far-fleld  detector  Is  d/N.  The 
achieved  resolution  will  be  XNF^/s.  Setting  these  two  expressions  equal 
and  using  the  fact  that  s  ■  2Nbs  gives 
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MFj  ■  2dA$/X.  (8-1) 

The  detector  elenent  spacing.  As,  therefore  has  an  leportant  effect  on 
several  other  paraneters.  Note  that  the  focal  length,  can  be 
Increased  to  conpensate  for  a  larger  d  and  that  Increases  In  the 
magnification  factor,  M,  can  be  used  to  decrease  F^. 

The  aperture  diameter,  0^,  of  lens  Lj  must  be  sufficiently  large  so 
that  rays  leaving  the  extreaws  of  the  target  and  traveling  In  directions 
corresponding  to  the  maximum  allowed  spatial  frequency  are  not 
vignetted.  The  maximum  spatial  frequency  Is  s/2XMFj  which  corresponds 
to  an  angle  6  with  sin  e  ■  s/2NFj.  The  aperture  must  be  of  diameter 
d  *  2Fjtan  e,  so,  for  e  small, 

Oj  •  d  ♦S/M.  (8-2) 

In  the  planned  experiment,  s  Is  less  than  d  and  N  Is  greater  than  1,  so 
the  aperture  diameter,  0^,  Is  slightly  larger  than  the  target  diameter, 
d. 


The  aperture  diameters  of  lenses  and  similarly  need  only  be 
large  enough  to  avoid  vignetting  and  can  be  easily  calculated.  The 
aperture  (see  Fig.  8>1)  In  the  common  focal  plane  of  lenses  and 
must  be  of  diameter  equal  to  the  far-fleld  detector  diameter,  s.  This 
ensures  that  the  Image  collected  by  the  Image  detector  has  the  same 
spatial  frequency  content  as  the  data  collected  by  the  far-fleld 
detector. 

The  magnification  from  the  target  to  the  Image  detector  Is  ^4/MF^ 
where  F^  Is  the  focal  length  of  lens  L^.  Since  at  least  2N  samples  of 
the  Image  must  be  detected,  F^  must  be  at  least 
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■  2KFjNAw/d 


(8-3) 


where  Aw  Is  the  detector  element  spacing  of  the  Image  detector.  By  an 
argument  similar  to  that  given  above  for  lens  L^,  the  aperture  diameter 
of  lens  must  be  at  least  s  >  w. 

One  possible  set  of  parameters  calculated  from  the  above  and  which 
represents  the  current  plan  for  the  active  experiment  Is: 

Target  diameter,  d  ■  25.7  mm 

Target  space-bandwidth  product,  N  -  128 

Lens  Lj  focal  length,  F^  ■  500  mm 

Lens  Lj  aperture  diameter,  0^  ■  27.0  mm 

Magnification  by  lenses  L^  and  L^,  M  *  6 

Far-fleld  detector  diameter,  s  ■  7.68  mm 

Far-fleld  detector  element  spacing.  As  >  30  microns 

Lens  L^  focal  length,  F^  ■  1000  nrni  (must  be  greater  than  896  mm) 

Target  Image  detector  diameter,  w  -  8.57  mm  (for  F^  >  1000  mm) 

Target  Image  detector  element  spacing,  aw  ■  30  microns 

Wavelength,  X  >  0.5145  microns. 

Light  Level  Consideration 

The  optical  Intensity,  I,  Incident  on  the  far-fleld  detector  will 
be  a  product  of  the  following  factors  (estimated  or  measured  values  are 
given  where  appropriate): 

Laser  power,  1  Watt 
Transmittance  of  spatial  filter,  0.5 
Transmittance  of  mask,  0.5 
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Reflectivity  of  target,  0.5 

Two-way  efficiency  of  beaiaspltter,  0.25 

Light  collection  efficiency  (assunlng  diffuse  target), 

Tr(s/2MFj)^/27r  - 

Transmittance  of  six  lenses,  0.9^ 

Gain  due  to  Image  Intensifler  (If  used),  6 

2 

Reciprocal  of  detector  area,  1/s 

The  resulting  expression  Is  approximately: 

1-2.1  G/(MFj)^  X  10“^  Matts.  (8-4) 

Using  the  values  from  the  previous  paragraph  gives  an  Intensity,  I,  of 
about  2.3  G  X  lO”  Hatts/cm  .  This  can  be  compared  with  a  laboratory 
measurement  with  a  system  similar  to  that  of  Fig.  8-1  (with  no 
Intensifler)  which  gave  an  Intensity  of  about  3  x  10  Hatts/cm  for  a 
cast  metal  target.  Since  current  solid  state  array  detectors  have  noise 

•ft  7 

equivalent  Intensities  of,  for  example,  about  1.4  x  10  Hatts/cm  for 
the  Fairchild  CC03000,  either  frame  Integration  or  an  Image  intensifler 
with  a  gain,  G,  of  at  least  100  Is  planned  to  be  used  In  the  active 
experiment. 

Software  Development 

Software  Is  currently  under  development  to  control  the  signal 
processing  hardware  and  to  permit  digitization  of  detector  signals, 
preprocessing  of  the  data,  and  computation  of  phase  retrieval  and  image 
reconstruction  algorithms.  An  outline  of  the  functions  planned  for  this 
experiment  control  program  Is  given  below. 
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1.  General 

Program  to  operate  In  command  file  and  Interactive  mode. 

Program  to  save  journal  of  all  commands  Issued  and  their  responses 
to  the  user  via  the  terminal  Including  comment  lines  In  order  to 
document  work  done. 

2.  Data  acquisition 

Digitize:  Digitize  and  store  image  In  Imaging  Technology.  Inc. 
(ITI)  frame  buffer  with  correction  for  calibrated  nonuniformities. 

Integrate:  Digitize  n  Images  and  sum  In  array  processor  (AP) 
with/without  normalization. 

3.  Image  transfer 

Transfer:  Move  Images  from  ITI  to  AP  and  hard  disk  and  between  AP 
and  hard  disk. 

4.  Image  display 

Display:  Display  AP  and  hard  disk  Images  on  ITI. 

Notes:  (1)  Conversion  from  32  bit  to  8  bit  data  needed 

(2)  Many  options  needed: 

Display  real.  Imaginary,  magnitude,  magnitude- squared, 
or  phase 

Apply  bias  and  scale  (as  In  y^ax-t-b) 

Display  absolute  value 

Magnify  by  2,4,8,...  (specify  subimage  to  be  displayed) 
Sample  to  give  256x256  Image  and  display  In  specified 
quadrant  of  ITI  display  (allows  four  Images  to  be 
displayed  simultaneously  for  comparison) 

Display  any  size  Image  In  any  location  of  display 

(3)  Some  of  these  options  can  be  done  by  altering  lookup 
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tables  In  ITI 

(4)  Values  above  and  beloM  the  8  bit  range  of  the  ITI 
should  be  clipped  at  0  and  255 

Llve/Ncmeory:  Toggle  between  displaying  video  Incoming  to  ITI  and 
data  In  ITI  franc  buffer. 

5.  Image  algebra  (all  In  AP) 

Add:  Add  two  Images. 

Subtract:  Subtract  two  Images. 

Multiply:  Multiply  two  Images. 

Divide:  Divide  two  Images  with  user  definable  result  for  divide  by 
zero. 

Scale:  Add  bias  and  scale  Image. 

.Threshold:  Hard  limit  above  and  below. 

Logic  operations  between  binary  Images. 

Magnitude:  Find  magnitude  or  magnitude-squared  of  an  Image. 

Phase:  Find  phase  of  an  Image. 

Convert:  Change  real  Image  to/from  complex  Image. 

Print:  Print  values  of  specified  small  part  of  an  Image. 

Statistics:  Find  mean,  variance  of  Image  and  magnitude-squared  of 
an  image. 

Maxmln:  Find  max  and  min  values  of  Image. 

Histogram:  Compute  histogram  of  Image  and  display  on  ITI. 

Convolve:  Convolve  Image  with  a  small  specified  convolution  kernel 
(allows  smoothing  and  other  operations  on  data). 

Interpolation:  Interpolate  from  one  sample  spacing  to  another. 

6.  Create  Images  (In  AP)  for  test  purposes 

Zero:  Zero  fill  an  Image. 

Create:  Place  a  rectangular,  circular,  or  triangular  region  of 


124 


specified  complex  value  at  a  specified  position  In  an  Image. 

Aperture:  Multiply  Image  by  a  binary  rectangular,  circular,  or 
triangular  aperture  located  at  a  specified  position. 

Noise:  Add  zero-mean  Gaussian  noise  with  specified  variance  to  an 
Image  (should  specify  seed  so  that  same  or  different  set  of 
random  numbers  can  be  generated)  (Also  Include  uniform  and 
Poisson  noise). 

7.  Image  warp 

Measure  warp:  By  use  of  calibrated  test  patterns,  measure  Image 
magnification  and  distortion. 

Remap:  Resample  Image  to  compensate  for  magnification  and 
distortion. 

8.  Iterative  algorithm 

Setup:  Allocate  and  load  Image  domain,  Fourier  domain.  Image 
domain  constraint,  Fourier  magnitude  constraint,  and  buffer 
arrays  In  AP. 

Iterate:  Iterate  n  times  using  specified  form  of  Iterative 
algorithm,  computing  and  printing  error  measures. 

Display:  Display  Intermediate  results. 

Save:  Save  results. 

9.  Image  error  computation 

Error  measure:  Compute  normalized  root-mean-squared  error  of 
complex  Image  or  of  Image  magnitude  relative  to  reference 
object,  taking  Into  account  Intensity  scaling  and  (for  complex 
Images)  constant  phase  shift. 
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10.  Help  Information 

Help:  Print  list  of  available  coamands. 

Print  help  Information  about  specified  coamiand. 

11.  Termination 

Stop:  Complete  all  commands  Issued  Including  writing  all  buffers 
to  hard  disk,  save  Journal  file,  and  return  to  UNIX. 

Journal  file  should  also  be  saved  at  each  step  In  case  of  program 
or  system  crash. 

8.2  PASSIVE  EXPERIMENT 

The  objective  of  the  passive  experiment  Is  to  demonstrate  Imaging 
(In  the  visible  or  Infrared)  of  a  noncoherently  Illuminated  or  emitting 
target  from  Intensity-only  or  Intensity  and  reduced  tolerance  phase 
measurements.  Several  candidate  experiments,  currently  under 
consideration,  are  described  below. 

Stellar  Speckle  Inteferometry 

Images  of  space  objects  from  ground  facilities  are  degraded  by  the 
effects  of  the  turbulent  atmosphere.  One  solution  to  this  problem  Is  to 
use  adaptive  optics  and  to  correct  for  the  phase  distortions  of  the 
atmosphere  In  real  time.  This  solution,  of  course,  requires  precise 
phase  measurement  and  compensation  In  real  time.  The  reduced  tolerance 
solution  Is  to  use  (stellar)  speckle  Interferometry  which  can  determine 
the  magnitude  of  the  Fourier  transform  of  the  target  Image  despite  the 
turbulent  atmosphere.  Phase  retrieval  and  Image  reconstruction 
algorithms  can  then  be  used  to  form  an  Image  of  the  target.  Specl  'ie 
Interferometry  has  been  used  by  astronomers  to  Image  simple  objects  such 
as  binary  stars.  The  technique  has  also  been  used  In  computer 
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siMulatlons  to  Image  more  complicated  objects.  An  experiment 
deannstratlng  the  use  of  speckle  Interferometry  to  Image  complicated 
objects  wuld  therefore  be  a  significant  step  forward  In  the  development 
of  reduced  tolerance  Imaging  through  turbulence.  The  experimental  setup 
required  could  use  much  of  the  equipment  from  the  active  experiment  with 
the  addition  of  glass  plates  with  small  thickness  variations  to  simulate 
the  effect  of  the  turbulent  atmosphere. 

Passive  Synthetic  Aperture  Imaging 

ERIM  Is  currently  developing  techniques  for  passive  synthetic 
aperture  Interferometric  Imaging  of  noncoherently  Illuminated  or 
emitting  targets  In  the  visible  and  Infrared.  These  techniques  require 
accurate  alignment  and  position  control  of  the  sensor  optics  (s1m11a>*  In 
degree  to  that  required  In  conventional  Imaging).  The  use  of  reduced 
tolerance  Imaging  techniques  could  reduce  these  accuracy  requirements. 
It  would  be  very  appropriate  and  effective  for  ERIM  to  Initiate  the 
research  to  combine  thes)*  two  techniques.  The  experimental  setup 
required  would  rely  heavily  on  the  equipment  used  for  passive 
Interferomtrc  Imaging.  However,  the  passive  Interferometric  Imaging 
program  plan  Is  such  that  Its  experiments  will  take  place  In  parallel 
with  or  after,  rather  than  before,  those  of  the  reduced  tolerance 
iMgIng  program,  so  It  may  not  be  feasible  to  link  the  two  experimental 
programs  together. 

Imaging  with  Phase  Diversity 

In  conventional  passive  Imaging  systems  where  the  Image  Is  degraded 
by  phase  aberrations  due  either  to  atmospheric  turbulence  or  to 
misalignment  of  segmented  optical  element  arrays.  It  Is  known  that  Image 
quality  can  be  Improved  by  using  Iterative  reconstruction  algorithms 
operating  on  two  2-D  Intensity  measurements.  In  the  case  of  turbulence, 
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these  tMO  meisurements  can  be  of  a  best  focus  Inage  and  an  Intentionally 
slightly  defocused  one  (that  Is.  a  quadratic  phase  error  Is 
Intentionally  Introduced  —  phase  diversity).  For  a  primary  mirror  made 
of  segments,  the  segments  may  be  slightly  moved  or  tilted  by  a  known 
amount  between  data  collections  In  the  Image  plane.  The  use  of  this 
approach  allows  reduced  tolerance  to  atmospheric  phase  or  to  accurate 
positioning  and  alignment  of  segmented  optics.  The  equipment  required 
here  would  again  be  similar  to  that  used  In  the  active  experiment  except 
that  piezo-electric  actuators  controlling  multiple  segments  would  be 
required. 

Further  planning  of  the  passive  experiment  Is  purposely  proceeding 
slowly  until  the  active  experiment  Is  approximately  one-half  complete. 
If  the  stellar  speckle  Interferometry  experiment  Is  chosen,  then  the 
optical  equipment  requirements  are  not  expected  to  Involve  any  major 
purchases  beyond  those  of  the  active  experiment  and.  In  any  case,  the 
electronic  signal  processing  hardware  already  purchased  will  serve  both 
experiments. 
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Appendix  A 

PROOF  OF  THE  UNIQUENESS  THEOREM 

In  this  appendix  the  uniqueness  theorem  presented  in  Section  3.2.3 
is  proven. 

Let  S,  T,  for  n  .  0,  .  .  S  -  1,  and  for 

n  ■  0,  .  .  T  -  1  be  defined  as  in  Section  3.2.4.  Therefore  all  the  q^ 

referred  to  in  this  appendix  are  reference  points.  Also,  let  U,  s^  and 

for  n  -  0,  .  .  .,  S  -  1  and  and  for  n  ■  0,  .  .  T  -  1  be 
n  n  •'n 

defined  as  in  Section  3.2.5.  Let  t^  be  the  side  of  [R(M)]  with  end- 
points  p„  and  P,„.,,^<|  j  (see  Figure  A-1 )  and  let  u„  .  P,„.,  T  '  "n 
for  n-0,  ...,T-1.  We  note  for  future  reference  that  for 
V,  we^^,  <v,  Uv>  -  0,  U^v«-v,  <Uv,  Uw>  •  <v,  w>,  and 
<v,  Uw  >  ■  <Uv,  U'^w^  ■  -  <  Uv,  w)  .  The  proof  of  the  uniqueness 
theorem  in  Section  3.2.3  requires  a  series  of  lemmas. 

Leflfma  A-1 ;  R(M)  is  a  mask  and  R(R(M))  .  R(M). 

Proof:  Suppose  it  can  be  shown  that  every  vertex  of  [R{M)]  is  opposite 
some  side  of  [R(M)].  Since  at  most  one  vertex  can  be  opposite  a  given 
side  and  the  nimber  of  vertices  equals  the  number  of  sides,  it  would 
then  follow  that  every  side  must  have  a  vertex  opposite  it  and  there¬ 
fore  no  two  sides  can  be  parallel,  hence  R(M)  is  a  mask.  Also,  since 
every  vertex  is  opposite  a  side,  R(R(M))  is  the  set  of  all  vertices  of 
[R(M)]  which  is  equal  to  the  set  R(M),  hence  R(R(M))  -  R(M).  Thus  it 
suffices  to  show  that  every  vertex  of  [R(M)]  is  opposite  some  side  of 
[R(M)]. 
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V4  =  P0  = 


FIGURE  A-1. 


=  P3  =  ^4 


THE  SET  [R(M)]  IS  THE  CONVEX  POLYGON  WITH  SIDES  t 
i  ■  0,  ....  4. 
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T  -  1.  Let  m  be 


The  vertices  of  [R(M)]  are  m  -  0,  .  . 
fixed  but  arbitrary,  0  £  m  £  T  -  1.  Then  is  also  a  vertex  of 
[M]  and  hence  p^^  -  Vj^  for  some  k.  Let  v^  be  the  vertex  of  [M] 
opposite  side  ^  of  [M]  and  let  v^^  be  the  vertex  of  [M] 

opposite  side  Sj^  of  [M].  Then  v^  and  v^^  are  in  R(M)  and  v^  -  p^ 
for  some  n.  (Refer  to  Figure  A-1  and  take  m  ■  0.  In  this  case  k  >  4, 
a  -  7,  b  -  1  and  n  -  2.)  4nd  v^  are  the  same  vertex,  then 

there  can  be  no  side  of  [M]  opposite  Vj^.  But  \  ■  P^  *  R(M) 
and  so  V|^  must  be  opposite  some  side  of  [M].  Therefore  v^  *  Vfa.  It 
then  follows  that  v^^  -  P(n+i)rBod  T*  ***  shown  that  p^^  is 

opposite  side  t^  of  [R(M)].  That  is,  we  wish  to  show  that  for 
0£j£T  -1  andj^m, 

<Pj,  Uu^  >  <  <p|j^,  Uu^  >  .  (A-1) 

Since  v^  is  opposite  side  5(|(_]),nod  S  follows  that  for 

0  <  i  <  S  -  1  and  i  ^  a. 


since  v^  is  opposite  Sj^ ,  for  0  £  i  £  S  -  1  and  i  4  b, 

<  UW|^,  v^  -  Vjj  >  <  0. 


(A-2) 


(A-3) 


Setting  i  -  b  in  (A-2)  we  obtain  <  )fl,od  S’  ^b  "  '^a  ^  ^  ^ 

thus 
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^'’(k-ljmod  S  ‘’m’  **“0^  “  ^ '(k-l jiwxl  s  '  *k’  *^^'’(n+1)mod  T  ”  •’n^^ 

■  -<''(k-1)m(xl  S*  “•'b  -  '»)> 

■  <‘*'(k-t)inod  S*  'b  -  'a> 


Setting  1 


<  0. 


a  in  (A-3)  we  obtain  <Uw.  ,  v 

K  A 


(A-4) 

V|^  >  <  0  and  thus 


^''(k^Diiiod  S  -  "«•  ““r,>  •  <’'(k«1)niO()  S  '  ’k-  ‘'<l>(n*1)nK)d  T  '  l>n>> 

■  U(y^  -  y,)> 

-  ''a  -  '^5  >  <  0. 

(A-5) 


since  5,  y^(  .  b„)  and  5  are  distinct  yertices 

of  CK).  the  yector,  3  .  .„p  ^  .  p^ 

linearly  Independent.  Let  p  .  R(1,).  p  *  3.  p^. 

Then  there  exist  real  ni«bers,  a  and  i,  such  that 


S* 


P  -  S  “  ®  ^''(lc♦1)tl»od  S  ■ 

Also,  since  <  ,  Uw^  >  <  <  p,  >  , 

0  <  <  P  -  ''ij. 

"  “  ^''(l(-l)(iiod  S  ■  '’m*  ^^''(k^Umod  S  " 

^  •^''(k^Dmod  S  "  Pm»  )iiiod  S  ‘ 

*  *^''(k-i)i"Od  S  ■  '’hi*  ^^''(k^l)mod  S  ‘ 
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(A-7) 


"  ■*^*'(k-l)mod  S’  S  "  ''k^^ 

-  5.  ^f^^4.])no6  S  ■  ''k^  • 

Since  5.  S  -  >  >  0.  it  follows  from  (A-7) 

that  a  >  0.  Similarly, 

0  <  <p  -  y^.  S> 

-  <P-P„.  U{p^-  S)> 

“  “^''(k-l)mod  S  '  '’m*  "  ''(k-l)mod  S^^ 

*  •^''(k*l)fflod  S  ■  ^’m’  ‘  ''(k-l)mod  S^^ 

-  -«<Pm  -  5.  U(p^  -  j)) 

"■  •<''(k^1)mod  S  ■  '*m’  "  ''(k-l)mod  S^> 

•  •^''(k>^l)mod  S  ‘  ^m’  ‘  ''(k-l)mo<J  S^> 

•  •^''(k>l)mod  S  ■  ''k’  '^(k-l)mod  S  ^  .  (A-8) 

<''(k>l)mod  S  -  ''k’  ^(k-l)mod  S  >  " 
that  a  >  0.  Using  (A-4),  (A-5),  (A-6)  and  the  fact  that  a  >  0  and 

.  >  0,  «  h.,e,  for  p  *  5.  P,.  5 

<P  -  P„.  Uu„>  .  •<»,k.,),„0  s  -  V  ''“n> 

•  •^''(krDi.oP  S  -  '’m‘  'J“n> 

<  0.  (A-9) 

Inequality  (A-1)  now  follows  from  (A-4),  (A-5)  and  (A-9).  This 
completes  the  proof  of  Lenna  A-1. 
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Lewwa  A-2;  The  number,  T,  of  points  In  R(M)  is  odd  and  if 
K  -  (T  -  1 )/2  and  0  £  j  <  T  -  1  then  pj,  as  a  vertex  of  [R(M)],  is 
opposite  side  j  of  [R(M)]. 

Proof;  It  follows  from  Leflwa  A-1  that  every  side  of  [R(M)]  has 
exactly  one  vertex  of  [R(H)]  opposite  from  It  and  every  vertex, 

I.e.,  every  point  In  R(M),  Is  opposite  from  exactly  one  side.  Thus 
there  Is  a  positive  Integer  K  <  T  -  2  such  that  P|^^^  Is  opposite 
side  t^.  Then  p^  Is  opposite  side  t^^^^  and  p^  Is  opposite  side 

‘(K*2)«io<l  T-  9«"«r*ny.  pj  (s  opposiu  side  Set- 

ting  j  •  T  -  K  we  obtain  Py_^  Is  opposite  t^.  But  Pj^^^  Is  opposite 
t^.  Therefore  T-K«K^lorT«2K^l.  This  completes  the 
proof  of  Lenwia  A-2. 

Lenswa  A-3;  For  x  c  [M],  x  ^  Pj,  y,j-0,  ...,T-1, 

<qj.  yj  >  <  <x.  yj  >  <  <R(j>i)mod  T»  ^  * 

Proof:  It  suffices  to  show  that  the  Inequalities  hold  for  all  ver¬ 
tices  V  of  [M],  V  ^  Qj,  y  Let  J  be  fixed  but  arbitrary. 

For  convenience  let  m  ■  (kj  -  l)mod  S  and  n  y.  Then 

■  ""m  -  “"n- 

First  we  will  show  that  j.  **  *  vertex  of  [M],  is 

opposite  side  s_  of  [M].  Let  v.  be  the  vertex  of  [M]  opposite  side 
s^  of  [M]  and  let  v^^  be  opposite  y 
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R(K)  and  >  p|^  for  some  k.  (Refer  to  Figure  A-1  and  take  j  ■  4. 
Then  k.  •  1,  m  ■  0,  (j  1  )mod  T  ■  0,  n  ■  k  ■  4,  a  ■  4,  b  ■  6  and 

J 

k  •  0.)  By  the  argument  In  the  proof  of  Lemma  A-1,  j 

"(JK)..),!  T  >  ■  "j  ■’'(..♦Dmod  $)•  «  * 
opposite  side  tj^  of  [R(M)].  By  Lemma  A-2,  j  opposite 

'(J*1)K  nod  T-  •>>  l-MPa  A-),  t,  -  TOd  T 

k  .  (J*1)K  nod  T.  Thus  T  ’  nod  T  *  ^k  '  %• 

therefore  Rjj^ijnod  T*  **  *  vertex  of  [M],  Is  opposite  side  s^^  of 

[H]. 

By  a  similar  argument  It  can  be  shown  that  q^,  as  a  vertex  of 
[M],  Is  opposite  side  s^  of  [M]. 

"(J^Dnod  T  '>»”'**  ‘n>  *  "(j^Dnod  T- 

^‘'(J*l)n»d  f  “*m>  >  <'•  '*'ni>  '»•  <'>(j»l)nod  T  *  »•  %>  > 
Also,  since  v  c  [M],  <v  -  v^,  Uw^>  >0.  Therefore, 

^^(j*l)mod  T’  ^  •  ^‘’(j^l)mod  T  '  ''»  > 

•  <‘^(j>l)mod  T  - 

"  ^^(J^l)mod  T  "  ^m^  *  ^‘^(j^l)mod  T  ‘  ^n  ^ 

•  <'»(j^l)mod  T  -  %>  -  <%  -  ‘^n^ 

•  <‘^(j-1)mod  T  -  l^n> 

>  0. 

(A-10) 
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Since  Qj  is  opposite  s^^  and  v  ^  Qj,  >  <v,  Uw^  >  or 

<qj  -  V,  Uw^>  >  0.  Also,  since  v  «  [M],  <  v  -  3,  Uw^  >  >  0. 

Therefore 

<V.  yj>  -  <qj.  yj>  -  <*  -  qj.  y^> 

•  <’  -  '>J' 

-  *  <qj  -  V.  Uw^> 

■  -  ’kj-  'J*m>  *  -  ’• 

■  -  ’'(hoDckkI  S*  -  '’•  '*'n> 

>0.  (A-11) 

It  noM  follows  from  (A-10)  and  (A-11)  that 

<qj.  yj>  <  yj>  <  <‘’(j+i)mod  T»  completes  the 

proof  of  Lermia  A-3. 

In  the  remainder  of  this  appendix,  all  modulo  arithmetic  will  be 
mod  T.  For  convenience,  we  define  m®  n  -  (m  ♦  n)mod  T. 

The  next  lemma  asserts  that  Oj  and  Oj  ^  ^  have  unique  separation 
in  M. 

Lewna  A-4;  For  0  <  j  ^  T  -  1 ,  if  x^ ,  x^  e  M  and  x^  -  x^  ■  Qj  q  ^  ”  ^j’ 
then  x^  ■  Qj  ^  ^  and  x^  ■  ^j* 

Proof;  If  either  x^  ^  •  1  ^2  ^  ‘Ij  follows  from 

Lemma  A-3  that 
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<Xi  -  X2.  yj>  -  <x^.  yj>  -  <X2,  y^  > 

<  <«>j  ®  1*  >'j> 

which  contradicts  the  assumption  that  x^  -  X2  ■  Qj  ^  ^  '  ^J* 

Therefore  x^  -  ^  ^  and  X2  ■  qj.  This  completes  the  proof 

of  Lenina  A-4. 

2 

Let  g  and  h  be  complex-valued  functions  on  Z  and  let  g  *  h 
denote  the  convolution  of  g  and  h.  That  is,  if  f  -  g  *  h  then 

f(x)  -  g(u)h(x  -  u).  (A-13) 

ueZ^ 

We  define  .S(g)  ♦  vS(h)  >  1  x  ♦  y;  x  e  J{g)  and  y  c  iS(h)  (  .  The  fol¬ 
lowing  Lemma  is  fundamental. 

Lemma  A-5;  If  f  -  g  *  h  then  [vS(f)]  •  [vS(g)  ♦  iS(h)]. 

Proof:  It  follows  from  (A-13)  that  5(f )  9  sS{g)  ♦  vS(h)  hence 
[vS(f )]  9  [^(g)  *  'S(h)].  It  remains  to  show  that  [vS{g)  vS(h)]  c 
[vS{f)]. 

Let  X  be  a  vertex  of  [»S(g)  ♦*S(h)].  Then  there  exists  a  y  c  31^ 
such  that  for  x*  c  J(g)  ♦  J(h)  and  x'  ^  x, 

<  x',  y>  <  <x,  y  >  .  (A-14) 

Also,  since  x  is  a  vertex  of  [J(g)  ♦  »S(h)],  it  follows  that 
X  c  J(g)  *  J(h)  and  hence  there  exist  x^  c  iS(g)  and  X2  c  3{h) 
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such  that  X  -  >  X2.  We  will  show  that  this  decomposition  of  x  Is 
unique.  Suppose  x  -  x'^  x'2  with  x'^  c  iS(g)  and  x'2  c  vS(h).  Then 

y>  *  <X2.  y>  •  <x^  ♦  X2.  y> 

-  <x.  y> 

-  <xj  ♦  xj.  y> 

-  <xj,  y>  +  <x^.  y>  .  (A-15) 

Therefore,  either  <x'^,  y>  >  <x^,  y>  or  <x'2.  y>  >  <X2,  y> 
or  both.  Suppose  <x'p  y>  2  y^*  Let  x'  ■  x'^  ♦  X2.  Then 

x'  e  J(g)  *  S(h)  and 

<x'.  y>  .  <xj.  y>  ♦  <X2.  y> 

>  <x^.  y>  ♦  <X2.  y> 

■  <  X,  y  >  .  (A-16) 

therefore,  by  Inequality  (A-14),  x'  ■  x  which  Implies  that 

x'^  a  x^  and  hence  x'2  -  X2.  If  <x'2,  y>  2  *  similar 

argument  leads  to  the  same  conclusion.  Therefore  the  decomposition 

X  -  x^  X2  with  x.|  c  J(g)  and  X2  c  «5(h)  Is  unique.  Now  suppose  for 

2 

a  particular  Uq  c  Z  ,  g(ug)  h(x  -  Ug)  ^  0.  Then 
Ug  c  3(g),  X  -  Ug  e  3(h)  and  x  -  Ug  ♦  (x  -  Ug).  By  the  uniqueness 
of  the  decomposition  of  x  It  follows  that  Ug  >  x^  and  x  -  Ug  >  X2. 
Therefore,  f(x)  ■  g(x.|)  h(x2)  ^  0  and  x  c  3(f).  Since  x  was  an 
arbitrary  vertex  of  [3(g)  ♦  3(h)],  It  follows  that  all  the  vertices 
of  [3(g)  ♦3(h)],  are  In  3(f)  and  therefore  [3(g)  ♦3(h)]  c  [J(f)]. 
This  completes  the  proof  of  Lemma  A-5. 
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We  are  now  ready  to  prove  the  theorem. 


Proof  of  Theorem;  Since  r  ■  r^ ,  it  follows  from  the  results  of 

Bruck  and  Sodin  [3.3]  that  there  exist  functions  g  and  h  with  finite 

2 

supports  and  a  vector  d  e  Z  such  that  f  -  g  *  h  and  f^(x)  ■ 

if  p 

g  *  h^(x  -  d)  where  h^(x)  «  h  (-x)  for  x  c  Z  . 

We  have  R(M)  c  ^(f)  c  ^(g)  +  iS(h).  Therefore  there  exist 
ag,  .  .  .,  ay_^  e  iS(g)  and  bg,  .  .  .,  by_.j  e  iS(h)  such  that 

*  **]*  '^  “  °*  *  *  **  ^  ‘  ’  *  ^ 

Now  let  J  be  fixed  but  arbitrary  and  let  x  c  iS(g)«  x  ^  aj.  We 
will  show  that  <aj,  yj>  <  <x,  yj^*  Suppose  to  the  contrary  that 
<  X,  yj  >  <  <aj,  Yj  ^  •  Let  x'  -  x  ♦  bj.  Then,  using  Lemma  A-5, 

X'  c  ^(g)  ♦  vS(h)  c  [j(g)  ♦  ^(h)]  .  [J(f)]  c  [M].  Also, 

<x',  yj>  •  <x,  yj>  ♦  <bj,  yj  > 

<  <«j.  yj>  *  <bj.  yj> 

•  <qj.  Yj)  .  (A-IB) 

It  now  follows  from  Lenwa  A-3  that  x'  -  q.  which  implies  that 

V 

X  ■  aj  contradicting  the  assumption  that  x  ^  aj.  Therefore 
^  ^  ^  ^  ^  *  similar  argument  it  can  be  shown  that 

if  X  c  ,S(g)  and  X  ^  aj  ^  ^ ,  then  <x,  Yj  >  ^  0  ]  •  ^ 

i  f  X  t  g )  and  X  ^  a  j ,  aj  0  ^ .  then 

<aj,  yj  >  <  <x,  yj>  <  9  i*  Yj  >  •  {A-19) 

Also,  by  similar  arguments,  it  can  be  shown  that  if  x  c  3{h)  and 
X  ^  bj,  bj  Q  ^  then 
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(A-20) 


<bj.  yj>  <  <x.  yj>  <  (bj*,.  yj>. 

Let  i5(g)  -  vS(h)  -  |  x  -  y:  x  e  J(g)  and  y  c  J{h)  f  .  Since 
J(h^)  -  -*S{h),»S(g)  ♦vS(h^)  ■  xS(g)  -  J(h).  Also,  since 
f](x)  ■  g  *  h^(x  -  d),  it  follows  from  Lemma  A-5  that 

[.S(f,)]  -  [J(g)  -  J(h)]  ♦  d.  (A-21) 

We  will  show  that 

aj  -  bj  ^  ^  ♦  d  e  M.  (A-22) 

We  have  ®  1  *  ‘^(9)  “  *5(h).  Let  j  be  fixed  but  arbitrary  and 

let  X  e  »S(g)  -  5(h),  x  Oj  -  bj  ^  Then  there  exist  x^  c  5(g) 

and  X2  c  5(h)  such  that  x  •  x^  -  X2.  Since  ^  ^  •  1*  either 

x^  ^  aj  or  X2  ^  bj  Q  ^  or  both.  In  any  case  it  follows  from  (A-19) 
and  (A-20)  that 

<x,  yj>  -  <Xp  yj>  -  <X2,  > 

^  ®  1’ 

■  ^  a  j  —  bj  ^  ^  *  {A~23) 

Therefore,  *j  "  0  1  is  *  vertex  of  [5(g)  -  5(h)]  and  by  (A-21), 

aj  -  bj  ^  ^  ♦  d  is  a  vertex  of  [5(f^)].  Therefore 
aj  -  bj  Q  1  +  d  e  5(fi )  9  M  and  (A-22)  fol1ov«s. 

By  a  similar  argument  it  can  be  shown  that 


Now 


aj  0  1  -  bj  >  d  .  M. 


(A-24) 
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(»j  -  Kj*,  *  d)  -  (dja  ,  -  bj  ♦  d)  .  (»j  *  bj)  -  (dja,  *  bja,) 

-dj-Qj#,.  (A-25) 

By  Lenina  A-4,  ®  1  *  ^  which  we  obtain 

bj  ♦  hj  ^  ^  ■  d.  (A-26) 

From  bj  ♦  bj  Q  ^  ■  d  and  bj  ^  ^  ♦  b j  q  2  ■  obtain  bj  ■  bj  ^  2* 
Since  by  Lenina  A-2  T  is  odd,  it  now  follows  that  b^  >  b^  -  .  .  ■ 

^T-1  using  (A-26),  we  obtain 

bg  ■  b^  ■  .  .  .  a  by  ^  ■  d/2.  (A-27) 

From  (A-20)  and  (A-27)  we  obtain  ,S(h)  a  { d/2  [  .  Therefore,  for  x  c 

f(x)  a  g  *  h  (x) 

a  h(d/2)g(x  -  d/2).  (A-28) 

If  h(d/2)  a  0  then  f  would  be  identically  zero,  contradicting  the 
assianption  that  R(M)  c  J(f).  Therefore,  h(d/2)  ^  0.  Now,  for 
X  e  Z^, 

f^(x)  a  g  *  h^(x  -  d) 

a  h*(d/2)g(x  -  d  ♦d/Z) 
a  h*(d/2)g(x  -  d/2) 

-  <»f(x).  (A-29) 

where 
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Since  la 


a 


it 


h  (d/2) 

fnw- 


■  1,  this  completes  the  proof  of  the  theorem. 
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Appendix  B 

PROOF  OF  PROGRAM  FOR  IMPLEMENTING  RECONSTRUCTION  ALGORITHM 

It  will  be  shown  in  this  appendix  that  the  program  presented  at 
the  end  of  Section  3.2.4  computes  f{d„)  for  n  -  T,  .  .  . ,  N  -  1 .  It 
will  be  assumed  that  f(R„)  for  n  -  0,  .  .  .,  T  -  1  has  been  computed 
as  described  in  Section  3.2.4.  Then,  since  0  £  m^  £  T  -  1,  f{n'p) 
has  been  computed  for  n  ■  T,  .  .  N  -  1. 

For  T  £  n  £  N  -  1  we  have 

r((ln  -  <l„  )  -L  f(y)f*(y  -  q„  *  q„  ).  (B-1) 

"  y.z2  " 

If  y  e  S{f)  and  y  -  c  J(f)  then,  since  iS(f)  9  M,  it  follows 

n 

that  y  c  M  and  y  -  q^  ^  q^^  c  M  or,  equivalently,  y  e  M  +  q^  -  q^^  . 

n  n 

Therefore, 

y  c  M  n (M  +  -  q^  )  C  I  .  (B_2) 

n 

Hence, 

n 

'■(‘In  -  ‘'m  )  -  E  -  ''n  *  "m  ) 

"  k.O 

n-1 

•  >  *  E  *  Om  >•  (B-3) 

"  k.O  " 

Thus, 
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H(n):  is  computed  correctly  for  0  £  j  £  n.  (B-5) 

By  the  derivation  in  Section  3.2.4,  H(T  -  1)  is  true.  Now  assume 

T  £  n  £  N  -  1  and  H(n  -  1)  is  true.  We  want  to  show  that  H(n)  is 

true.  It  suffices  to  show  that  f(Q„)  is  computed  correctly.  Let 

all  variables  have  the  values  that  they  have  at  Step  2  of  the  pass 

through  the  loop  in  which  f(q^)  is  computed.  It  must  be  shown  that 

all  values  of  f  appearing  in  the  right-hand  side  of  Eq.  {B-4)  are 

correct.  Since,  by  assumption,  H(n  -  is  true  it  follows  that 

f(q(^).  k  -  0,  .  .  .,  n  -  1,  have  the  correct  values.  Also,  as 

mentioned  above,  f(qj^  )  has  the  correct  value.  Now  let 

n 

X  ■  with  0  £  k  £  n  -  1.  If  x  ^  M,  then  f(x)  ■  0.  In 

n 

this  case,  since  v5{f)  9  M,  it  follows  that  x  f  ^^(f)  and  therefore  0 

is  the  correct  value  for  f(x).  Now  assume  x  e  M.  We  have 

^  therefore  *  *  ^  ~  q^,  •  Hence 

n  n 

X  c  Mn(M  -  q^  q^  )  9  I  qQ,  .  .  .,  qp_^  }  ^  follows  from  the 

n 

induction  hypothesis,  H(n  -  1),  that  f(x)  has  the  correct  value  when 
the  value  for  f(q„)  is  computed.  Therefore  f(qp)  is  computed 
correctly  and  H(n)  is  true.  By  induction  it  follows  that  H(N  -  1) 
is  true  and  hence  f(n)  is  computed  correctly  for  n  -  0,.  .  .,  N  -  1. 
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Appendix  C 

PROOF  OF  THE  ALGORITHM  FOR  GENERATING  RECONSTRUCTION  ALGORITHMS 

In  this  appendix  it  will  be  shown  that  the  program  in  Section  3.2.5 
generates  a  reconstruction  algorithm.  First,  it  will  be  shown  that 
the  loop  is  not  infinite  and  hence  the  program  produces  sequences 
dy,  .  .  and  nij.,  .  .  m^_.j.  Secondly,  it  will  be  shown  that 

IF  q  ■  (Pq,  .  .  .,  qfj_^)  and  m  •  (tiiy,  .  .  .,  then  (q,  m)  is  a 

reconstruction  algorithm. 

In  the  following,  0  will  be  used  to  denote  both  the  number  zero 
2 

and  the  origin  of  31  .  Context  should  prevent  any  confusion. 

If  X,  y,  z  let  [x,  y,  z]  denote  their  convex  hull  in  51^. 

If  X,  y  and  z  are  non-col  1  inear  then  the  interior  of  [x,  y,  z]  is 
given  by 

int[x,  y,  z]  ■  lax  *  by  *  cz:  a,  b,  c  >  0  and  a  ♦  b  ♦  c  -  1  1  . 

(C-1) 

Then  0  c  int[x,  y,  z]  if  and  only  if  x,  y  and  z  are  non-col  1  inear 
and  there  exist  strictly  positive  numbers  a,  b,  c  such  that 
ax  by  *  cz  -  0.  (The  sum  of  a,  b  and  c  can  be  made  equal  to  1  by 
dividing  each  of  these  numbers  by  their  sum  if  necessary.) 

The  following  lemma  will  be  needed. 

2 

Lemma  C-1;  If  Up,  e  51  ,  n  ■  1 ,  2,  3,  Vp  ^  >0 

<Mp,  >  <  0  for  n  ^  m,  then  0  e  intCu^,  »3] 

0  c  int  [vp  V2,  vj]. 
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Proof:  By  s>»iinetry  it  suffices  to  show  that  0  c  int[u^,  U2i  u^]. 


First,  we  will  show  that  up  U2  W3  are  non-col  1  inear.  If 

they  were  coll  Inear  then  one  of  them  would  be  in  the  convex  hull  of 

the  other  two.  Say  is  in  the  convex  hull  of  W2  >*3* 

there  would  exist  numbers  T2  and  Tj  such  that  T2,  ^  Ot  "^2  "^3  ■  ^ 

and  ui  ■  12^2  *  ^3“3*  v.j>  •  ^T2M2  '^3^31  ■ 

T2<m2i  *^3  ^^3*  ''l  ^  ^  contradicting  the  assumption  that 

<Mp  v^>  >0.  Therefore  u^,  U2  and  are  non-col  linear. 

2 

Since  any  three  vectors  in^  are  linearly  dependent,  there 
exist  three  numbers  o^,  02  and  03,  not  all  zero,  such  that 

♦  a2U2  ■*“  ®3W3  -  0.  {C-2) 

Since  < u„,  >  >  0,  u„  ^  0,  n  •  1 ,  2,  3.  Therefore  no  two  of  the 

n  n  n 

o„  can  be  zero, 
n 

Now  suppose  -  0.  Then  02  ^  0  ^  03  and  it  follows  from  (C-2) 

that 


hence 


0  <  <  U2.  '»2  ^ 


Since  <113,  V2>  <0,  it  follows  from  (C-4)  that 


{C-4) 
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{C-5) 


Also, 


!i 

"2 


< 


0. 


0  >  <U2.  ^  ^**3*  ’  (C-6) 

and  since  (uj,  >  <0,  it  follows  from  (C-6)  that  -03/02  >  0  which 

contradicts  (C-5).  Therefore  ^  0  and  by  symmetry,  02  ^  0  ^  03. 

By  multiplying  the  o^'s  by  -1  if  necessary,  we  may  assume  that  >  0. 
Now 

«1  <u^,  >»i>  *  <m2.  '»!>"'  ®3  <‘*3*  '*!> 

.  <a^ui  ♦  «2‘‘2  ^  V3»  ''!> 

.  0.  {C-7) 

Since  >  0  and  <u^,  v^>  >  0, 

»2  <  ‘‘2»  ®3  <“3»  ''1  ^  ■  -O]  <  W].  '']  >  <0*  (C-8) 

Since  <  U2»  >  <0  >  <0,  at  least  one  of  the  numbers 

02  and  03  must  be  strictly  positive.  By  symmetry,  we  may  assume 
without  loss  of  generality  that  <^2  >  0* 

®1  <‘‘1*  '’3^  ®2<‘‘2»  ''3>  *  ®3<‘'3’  ''3> 

-  ^  a2‘‘2  ®3‘'3*  ''3  > 

-  0.  (C-9) 

Since  o-j  >  0,  02  >  0,  <up  v^)  <0  and  <m2.  '>3^  <  0»  IT  follows 
that 
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®3<‘'3’  ''3>  -  ''3>  ■  ®2<‘'2*  ''3> 

>  0.  (C-10) 

Since  <m2»  ''3^  >0*  follows  that  >  0.  By  the  coeinent  pre- 
ceedlng  the  letima,  it  now  follows  that  0  c  int[M^,  w^*  >*3^* 
completes  the  proof  of  Lenma  C-1. 

One  more  lenma  is  needed  before  proving  that  the  loop  is  not 
infinite.  As  in  Appendix  A,  we  define  j  ®  k  -  (j  ♦  k)mod  T. 

Lemma  C-2:  For  j  -  0,  .  .  T  -  1  and  k  -  2,  .  .  T  -  1, 

0  c  intCyj,  yj  ^  p  (-1) 

Proof;  The  proof  will  be  by  induction  on  k.  Let  k  -  2.  We  want  to 
show  that  0  c  int[yj,  Yj  (j  yj  9  2]*  Let 

“1  •  ®  1  ‘  ®  3’  “2  "  ‘’j  »  2  “  ®  r  “3  "  ®  2’ 

''1  •  ''2  •  -^j®  !•  '’3  •  2  • 

(C-11) 

By  Lemma  A-3,  <u„,  v„  >  >  0  and  <u„,  v„  >  <0  for  n  ^  m.  There- 
•'  ’n  n  n  m 

fore,  by  Lemma  C-1,  0  c  1nt[v^,  V2,  v^]  -  int[yj,  yj  q  1.  yj  ®  3]- 

Now  let  3  £  k  £  T  -  1  and  assume  the  lemma  is  true  for  k  -  1.  We 

want  to  show  that  0  c  int[yj,  ®  k^* 

that  0  c  int[yj,  y j  ®  ] »  yj  9  2^  therefore  there  exist  strictly 
positive  numbers  a-j,  02  03  such  that 

®l^j  *  ®2^j  ®  1  Vj  ®  2  •  °* 
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Applying  the  lenina  for  k  -  1  with  j  replaced  by  j  ®  1 ,  if<e  have 

k  1 

0  c  int[yj  9  p  yj  ®  2*  9  k^  therefore  there  exist 

strictly  positive  numbers,  t2,  Tj  such  that 

9  1  *  9  2*^  9  k  ■  (C-13) 

Multiplying  (C-13)  by  O3/T2  subtracting  the  result  from  (C-12) 
we  obtain 

^■jyj  *  ^2^j  9  1  ^  ^3  9  k  "  (C-14) 

whcpc  ■  (J ^  I  x^  •  ^2  *  ^3^3^ ^2*  x^  y  0 

and  x^  >  0.  We  will  show  that  X2  >  0.  From  (C-14)  we  obtain 

^2^j  9  1  *  ”^l^j  *  ^3^”^^  9  k*  (C-15) 

We  consider  two  cases. 

Case  1;  k  is  odd.  Then  k  -  1  is  even  and  by  (C-15) 

^2'^j  9  1  *  ”^l^j  *  ^3'^j  9  k’  (C-16) 

and,  using  Lemma  A-3, 

^2  9  k  9  1  ■  ‘'j  9  1’  9  1  > 

-  -Xi  <  <1j  ®  k  ®  1  •  ‘Ij  9  1  ’  ^  *  ^3  ^  *^3  9  k  9  1  ■  ‘’j  9  T  9  k  ^ 

-  <  q j  ©  1  "  9  k  9  1  *  ^  *  ^3  ^  9  k  9  1  "  ‘’j  9  1  ’  9  k  ^ 

>  0. 

(C-17) 
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Since,  by  Lenma  A-3,  ^  Q  |(  3  ]  ”  Qj  ®  i »  Yj  9  ]  >  >  0,  it  fol  lows 
from  (C-17)  that  >  0. 

Case  2:  k  is  even.  Then  k  -  1  is  odd  and  by  (C-15) 

®  1  "  '^l^J  ■  ®  k*  (C-18) 

and,  using  Lemma  A-3, 

^2  •  k  ■  •  1* 

-  x^  <qj  ^  1  -  Qj  3  ,5.  yj  >  ^  3  1  ‘  ‘’j  9  k’  9  k  ^ 

>  0. 

(C-19) 

Since  by  Lemma  A-3,  ^  9j  3  -  dj  3  yj  3  ^  >  >  0,  it  follows  from 

(C-19)  that  x^  >  0. 

U 

It  remains  to  show  that  y^,  yj  3  ^  and  (-1)  yj  3  |j  are  non- 

collinear.  Since  x^  X2,  x^  >  0,  It  follows  from  (C-14)  that 

0  e  [yj,  yj  3  (-l)Sj  3  ij].  Therefore  if  yj,  yj  3  ^  and  (-l)Sj  3 

are  col  linear  then  they  must  all  lie  on  a  line  through  the  origin. 

But  since  we  have  already  shown  that  0  c  int[yj,  yj  3  yj  3  j]* 

yj  and  yj  3  ^  cannot  lie  on  a  line  through  the  origin.  Therefore 

1/ 

yj,  yj  3  ]  and  (-1)  yj  3  are  non-coil  inear  and  hence 

1/ 

0  c  int[yj,  yj  3  i»  (-1)  yj  3  1^]*  T^»is  completes  that  proof  of  Lemma  C-2. 

In  order  to  prove  that  the  loop  is  not  infinite  it  will  be  shown 
that  the  parameter  n  in  the  program  in  Section  3.2.5  can  fail  to  be 
incremented  on  at  most  T  -  2  consecutive  passes  through  the  loop. 

The  proof  will  be  by  contradiction.  Accordingly,  assume  that  n  is 
not  incremented  on  T  -  1  consecutive  passes  through  the  loop. 
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Let  k  and  #  have  the  values  that  they  have  at  Step  2  of  the  first 
of  these  T  -  1  passes.  Let 

b^«m1nij:  O^j^N-T-l  and  j  ^  ^  ^  (C-20) 

for  a  •  0,  .  .  T  -  2.  Then 


- ' 

A 


(C-21) 


and 


•  a  ■  \  •  a.b^^  "  ^ 
for  a  ■  0,  .  .  T  -  2.  Let 

®  a.b  •  a^^^k  •  a»b  ^  ® 


(C-22) 


k  ®  a  ■  ^k  •  a,b, 


otherwise. 


(C-23) 


Then 


\  •  a^*a^  "  l  \  ®  a^'^k  •  a,b^ 


(C-24) 


If  X  c  and  4(x)  -  1,  then  x  c  D  and  x  -  dj^  ^  ^  j  for  some  j  2  ^a’ 

Therefore  Ih,,  *,(x)|  .  1^  «  ,(d|<  ,  j)l  <  ,  ("k  •  e.b  >1  * 

2  * 

\  9  a^*a^'  '*'^“**  ‘  Z  » 

«(x)  -  1  — >  lh,^Q^(x)l  <  \  ®a(Xa)-  ^^-25) 

Claim;  For  a  ■  0,  .  .  .,  T  -  2, 


(-l)*<x^  -  Qk  ^  ,,  yk  0  {T_i)>  >  0- 


(C-26) 
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Proof  of  Claim;  First,  we  will  prove  the  claim  for  a  -  0.  Since  by 


(C-21),  {C-22)  and  (C-23),  -  Xq)  .  1,  it  follows  that 

-  Xq  c  0.  Therefore,  by  Lemma  A-3, 

<\  -  *0*  J'k  •  (T-1)^  ^  ^k  •  (T-1)^  ’ 

or,  since  ^  -  ‘’k  ^  \  ®  1  • 

^*0  “  ^k  •  r  ^k  •  (T-1)^  ^  °* 

Thus,  the  claim  is  true  for  a  -  0.  Now  let  1  £  a  £  T  -  2  and 
the  claim  is  true  for  a  -  1.  By  (C-21),  (C-22),  and  (C-23), 
^^*k  ©a  ~  *a^  “  ^  hence  using  (C-25), 

\  •  (a-l)K  •  a  ■  *a^  i  K  •  (a-l)^“k  •  a  “  *a^l 

-  \  •  (a-l)^*a-l^’ 

or  equivalently, 

^*a-l  “  *k  «  a  *  *a’  ^k  e  (a-1 )  ^  ^  °- 
Similarly,  ^(Xj_^)  ■  1  and 

■’’k  •  I  ''k  ®  a^'a-l  ^  I 

i  \  •  a'S’ 

■  ®  a^®k  •  a  "  *a^’ 

and  therefore, 

^k  ®  a^*k  ®  a  ■  *a^  -  \  •  a^*a-l  ^  ’ 

or  equivalently. 


(C-27) 

(C-28) 

assume 


(C-29) 


(C-30) 


(C-31) 

(C-32) 


152 


"  ®k  ®  a  *  *a’  ®  a  ^ 

By  Lemma  C-2,  0  c  int[y,^  ^  y,^  ^  (-1)^"*  y,^  9  (j_ 

therefore  there  exist  strictly  positive  real  numbers  o-j,  o 
that 

®l^k  ®  (a-1)  *  ®2^k  ®  a  ^  *  ^k  ®  (T-1)  “ 

Since,  by  Lemma  A-2,  T  is  odd,  (-1)^“*  ■  -{-!)*,  and  from 
obtain, 

03(-1)*  ^k  ®  (T-1)  "  ®l^k  ®  (a-l)  "  Vk  ®  a* 
By  {C-35),  {C-30)  and  (C-33), 

’  “k  ®  a  *a*  ^k  ®  (T-1)^ 

■  ®1  <’‘a-l  "  ®k  ®  a  ^a’  ^k  ®  (a-1)^ 

"'^2<Vl  “  “k®  a  ^  J'k®a> 

>0, 

and  since  >  0, 

^*a-l  "  ®k  ®  a  ^  *a’  ^k  9  (T-1)^  - 
Substituting  %  ©  ^  •  <1|(  0  ^  \  ^  9  1  and  using  (C-37) 

induction  hypothesis. 


(C-33) 

1)]  and 
2 ,  such 

0.  (C-34) 

(C-34)  we 

(C-35) 


(C-36) 

(C-37) 
and  the 
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{-1)^<X^  -  q,^  0  ^  0  9 

>  -{-!)*  <  ®  a*  \  3  (T-1)^ 

-  (-1)^  ’  <Xj_^  ■  ®  a*  ^k  e  (T-1)^ 

>  0. 

(C-38) 

This  completes  the  proof  of  the  claim. 

Now  set  a  -  T  -  2  in  (C-26).  Since,  by  Lemma  A-2,  T  -  2  is  odd 
we  obtain, 

<XT-2'  ^k  ®  (T-1)^  <  ^<^k  •  (T-l)’  ^k  •  (T-1)^  • 

It  now  follows  from  Lemma  A-3  that  Xj_2  i  D  and  therefore 
^(xy_2)  ■  0  which  contradicts  either  {C-21)  or  (C-22).  Therefore  n 
cannot  fail  to  be  incremented  on  each  of  T  -  1  consecutive  passes 
through  the  loop.  This  completes  the  proof  that  the  loop  is  not 
infinite. 

It  now  follows  that  the  program  produces  sequences 
qy,  .  .  q|^_^  and  my.  .  .  .,  It  remains  to  show  that  if 

q  ■  (qQt  •  •  ••  d  m  ■  ( my .....  m|^  ^ )  then  (q.  m)  is  a 

reconstruction  algorithm  for  the  mask  M. 

We  have  R(M)  -  |  qg,  .  .  .,  qy_^  [  .  Let  T  n  <  N  -  1  and  let 
all  variables  have  the  values  that  they  have  after  Step  5  and  before 
Step  6  of  the  pass  through  the  loop  in  which  q^  and  m^  are  de¬ 
fined.  By  the  definition  of  b  in  Step  1  of  the  loop,  for  x  c 
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(C-40) 


If  X  e  0,  then  before  entering  the  loop  i  had  to  have  the  value  1  at 
X.  Thus  if  the  current  value  is  4{x)  -  0,  ^  must  have  acquired  the 
value  0  at  X  on  some  preceed^ng  pass.  That  is,  x  ■  q^,  for  some 
n'  <  n.  Therefore  x  e  |  q^,  .  .  q^_^  |  .  Thus,  for  x  c  D, 

^(x)  ■  0  xe  jqQ,  ...,  qp_^  I  •  (C-41) 

First,  it  will  be  shown  that 

Mn(M  ♦  q„  -  q^  )  C  {q^ . q J  .  (C-42) 

n 

Let  X  c  Mn(M  ♦  q^  -  q^  ).  We  want  to  show  that  x  e  |  qQ,  .  .  .  ,q^  [  . 

n 

If  X  e  R(M)  then,  since  n  ^  T,  we  are  done.  Since  q^^  -  dj^  1^.  if 
X  ■  dj^  ^  we  are  done.  Now  assume  x  i  R(M)  and  x  ^  dj^  1^.  Since 
X  ^  R(M),  X  e  0. 

Claim; 

<  l\(x)l  .  (C-43) 

Proof  of  Claim;  The  proof  of  the  claim  will  be  divided  into  two 
cases. 

Case  1; 

h,^{d^,b)  >  0.  (C-44) 

Since  %  ■  dj^  ^j,  it  follows  from  Step  5  of  the  loop  that 
m^  -  k.  Let  z  ■  x  -  dj^  ♦  qj^ .  Since  *  *  ^  ■*■  \  -  qj^,  it  follows 

that  z  c  M.  Also,  since  x  i  z  ^  Therefore,  by  Lemma  A-3, 


155 


<qi<,  yi^>  <  <2.  y|<> 

•  <X.  y,,)  -  y,,)  *  <q|,.  y,,)  .  (C-45) 

or  <d|^  y|^>  <  <x,  yj^  >  .  It  follows  that  h|^{d|^  <  h|^(x)  and 

since  h|^(d|^  >  0,  the  claim  follows. 

Case  2: 

h^(d^^b)  <  0.  (C-46) 

In  this  case  m^  «  k  ®  1.  Let  z  -  x  -  dj^  b  ‘Ik  9  1  •  Since 
X  c  M  +  d|^  b  ■  ‘Ik  9  1*  follows  that  z  e  M.  Also,  since  *  ^  d|^  b» 
z  ^  ^  ^ .  Therefore,  by  Lemma  A-3, 

<’'•  J'k>  -  <‘‘k,b>  A>  *  <''k«i-  ^k> 

-  <'•  j'k> 

<  <‘Ik  9  1»  ^k^  ’ 

or  <x,  y|^>  <  <d|^  b»  yk^  '  follows  that  h|^(x)  <  h|^(d|^  b)  ^‘’‘1 
since  lj)  <  0,  the  claim  follows.  This  completes  the  proof  of 

the  claim. 

Now  it  follows  by  implication  (C-40)  that  0(x)  «  0  and  since 
X  e  D,  it  follows  from  implication  (C-41)  that 
X  c  {qq,  .  .  .,  I  c  jq^,  .  .  .,  q„l  .  This  completes  the 
proof  of  (C-42). 

It  remains  to  show  that 

Mn(M  -  ‘Ik  b  ^  ‘^m  ^  •  •  ••  Vl  ^ 

*  n 
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Let  X  c  M  n  (M  -  d|^  ).  We  want  to  show  that 

’  n 

X  €  I  q^,  .  .  [  .  Since  n  ^  T,  if  x  c  R{M)  we  are  done.  Now 

assume  x  i  R(M).  Then  x  c  0.  The  proof  will  be  divided  into  two 
cases. 

Case  1 :  h|^(d|^  >0. 

In  this  case  -  k.  Let  ^  5  *  R|(*  Since 

X  c  M  -  d|^  b  *  follows  that  z  e  M. 

If  z  .  q,,  ,  then  X  .  »  q,,  *  ,  -  No»  if 

^(x)  >  1,  then  by  Step  2  of  the  loop,  q^  would  not  have  been  de¬ 
fined  on  this  pass  contrary  to  the  assumption  that  it  was.  There¬ 
fore,  ^(x)  -  0  and  by  implication  (C-41),  x  e  |  qQ,  .  .  .,  q^_^  |  . 

Now  assume  z  ^  qj^  ^  ^ .  Then  by  Lemma  A-3, 

<’'■  *  ^'“k.b-  ■I'k^  -  ^'"k-  J'k^ 

■  <*•  )'k> 

^  ^  ^k  ®  1’  '^k  ^ 

Recalling  that  •  (q|^  ^k  ®  \  follows  from  (C-49)  that 

^‘’k.b  -  “k-  J'k^  ^  *  ‘k’  )'k> 

or  h|^(d|^  jj)  <  -h|j(x)  and  since  h|^(dj^  >  0,  it  follows  that 

ih|^(dk  jj)|<  |h|^(x)|.  Now  by  implication  {C-40),  i(x)  ■  0  and  by 
implication  (C-41),  x  e  |  qg,  .  .  .,  q„_^  }  • 
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Since 


In  this  case  -  k  ®  1.  Let  z  ■  x  ♦  dj^  0  ^ . 

xeM-d|^jj-*-q|^01,  it  follows  that  z  e  M. 

If  z  -  q,^  then  x  -  qj^  ♦  qj^  0  ^  b  *  ®k  ‘  b 
argument  given  in  Case  1  above,  4(x)  «  0,  and  by  implication  (C-41), 
X  c  I  qQ,  .  .  .,  q^_^  I  . 

Now  assume  2  ^  qj^.  Then  by  Lemma  A-3, 

<q|j.  yi(>  <  <2.  y|j> 

-  <»•  ^k>  *  -  <''k»i-  j'k> 

from  which  it  follows  that  5)*  Since  jj)  <  0, 

it  follows  that  ^j)  I  <  lh|^(x)i.  Hence  by  implication  (C-40), 

^(x)  ■  0  and  by  implication  {C-41),  x  e  jqQ,  .  .  .,  q^_^  |  .  This 
establishes  (C-48)  and  completes  the  proof  that  (q,  m)  is  a  recon¬ 
struction  algorithm  for  the  mask  M. 


APPENDIX  D 


PARAMETER  ESTIMATION  AND  THE  CRAMER-RAO  LOWER  BOUND 

D.l.  INTRODUCTION 

The  eetimetion  probkai  we  eomsider  b  to  cetiBete  %  pirmmeter  vector 
j  »  (ii,  .  .  . ,  ei, )  froB  >  BeasnreBeat  vector  E.  ~  (A|»  •  •  •  >  KyY  The  estimate  of 
j  is  denoted  A  It  b  aesamed  that  the  conditional  probability  density 

p(l  Ii)  b  known.  We  consider  two  estimatiott  eavironBeats  that  differ  primarily  in 
their  prior  knowledge  of  j . 

A  word  OB  notation  b  in  order.  In  general,  we  will  ase  npper  ease  letters  to  denote 
random  variables  or  vectors  and  will  use  lower  case  letters  to  denote  their  possible  values 
and  also  most  noB*raadoB  variables  or  vectors.  Vectors  will  be  nnderscored.  Thus,  the 
received  vector,  which  contains  randoaness,  b  denoted  E-  A  function  of  a  random 
quantity  is  also  to  be  considered  random;  hence,  the  estimate  A  >9  capitalised.  As 
described  below,  the  parameter  vector  to  be  estimated  may  be  random  or  not.  In  the 
previous  paragraph,  we  elected  to  use  the  lower  case  g.  The  probability  density  of  a 
random  vector  such  as  E  *ill  he  denoted  Pgfg),  or  simply  p(i.),  when  no  confusion  can 
arise.  A  conditional  probability  density  such  as  that  of  E  given  A  will  be  denoted 
Pfi  l-i)  or  p{z,  |i).  The  average  of  a  random  quantity  will  be  denoted  either 
E  [A  I  or  [s  ].  The  latter  is  especially  helpful  for  expressions  like  E^  [dp^  (o)/dal 
which  otherwise  would  have  to  be  written  E[dpj^  {A)/dA\ 
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EDvlronmrat  It  NoarMdoin  ParanMtan 


Here,  the  parameter  rector  to  be  estimated  is  a  fixed  bat  aaknowo  rector  j ,  aod 
the  qoalitj  of  aa  estimator  is  jodfed  hj  tbs  mesa  sfnsrcd  crrsri  (MSE) 

-  •<)*!.  •  -I,i,  ...,L  .  (D-l) 

Thert  maj  bs  ao  ecrtifiably  best  estimator.  Tbs  nmsmmmh  likeftAeed  (ML)  estimator 
chooses  ifpCL|l)i*  largest  amoag  all  choicss  for  g .  Aa  estimator  is 

«ii4issed  if  E[A,]  ~  s,,  i  «  1,  .  . .  ,  L.  Tho  ML  estimator  might  or  might  not  be 
aa  biased. 

Earlroamaat  St  Raadon  Paramotors 

Here,  the  parameter  rector  to  be  estimated  is  a  raadom  rector  A.  vitb  knowa  pro* 
babilitj  deositj  p  (i ).  The  qaalitjr  of  aa  estimator  is  jodged  by  the  meaa  squared  errors 
(MSE) 

EU,  -A,f\.  I  -  I,  ....L  .  (D-2) 

The  best  estimator,  i.e.,  the  minimum  mean  sqaered  errer  (MMSE)  estimator  is 

i(£)-£ld  Ifi  -lI  (D*3) 

Ao  estimator  is  vniieecd  if  E(A,]  «  E(A,  ],  i  *  1,  .  .  . ,  L .  The  MMSE  estimator  is 
unbiased.  In  either  enrimnmsnt,  there  is  a  Cramer-Rao  lower  bound  (CRLB)  to  the 
MSE. 


Enrlronmoat  It  For  unbiased  estimator, 

-  s,)*l  >  ,  i  -  1,  .  .  .  ,L  (0*4) 
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(D.6») 


—  ^a| ^  f  U  1 1)  1“  f(i.  I  j ) 

where  the  expectations  aTerage  over  R, . 

Eaviroamont  St  For  estimators  hariag  a  etrtmn  “aniiMc^/ifte'’ property, 

. L  ,  (D-0) 

where 

(D-7) 

f<„  -  •■fU)  I^fU)!  (D-8a) 

In  this  eDTironment,  the  expectatioas  arerage  orer  both  R  and  ± .  la  the  scalar  case 
(i.e.,  L  ■*!)  the  "onbiased-like"  property  which  the  estimator  must  satisfy  is 

^n»^F4(«)|£'a(>^(i.)i  -  *1  —  0  •  (D-9) 

Notoot 

( 1 )  The  existence  of  the  shore  deriratires  is  presumed. 

(2)  The  matrices  J  and  K  are  functions  of  the  likelihood  distribution  p(£  |  ii )  and 
the  a  priori  distribution  p  (i ),  respectirely.  J  is  related  primarily  to  the  likelihood  dis¬ 
tribution  p(z.  I  A )  and  secondarily  to  the  a  priori  distribution  p(i ). 
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lo  the  CMC  where  ^  is  c  scslar  (L  «  1).  thcM  bounds  reduce  to  the  followinf. 


Euviroumont  li 


Environment  St 


EM  -.4)»1> 


d*  I 
— 

dc* 


-  . >«  I  ^  (r  I  • )  I  - 1^4 1  F  ( • ) 


Note  thnt  in  Environment  2  one  mnj  use  the  bounds  in  Eq.  (D«ll)  for  the  vector 
esse  M  well;  i.e., 

£M,  -\f\>  ([ «!; '’'•’f)) 

This  bound  is  presumsbly  simpler  to  compute  but  noi  m  good  m  thst  of  Eq.  (D>0). 


D.8  ADDITIVE  NOISE  PROBLEMS 


In  this  section,  we  focus  on  the  following  “sdditive  noise"  estimstion  problem  sod 
specisi  cMes  thereof: 

+  (D-13) 

where  /  is  n  known  function  ud  JH,  ^  yN^,  ,  Vy  is  s  random  noise  vector, 
independent  of  ^ ,  with  densitj  p^lu  )•  The  conditionsi  densitj  of  E  given  A  ~  1  is 
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I  - /U))  • 


(D.14) 


D.3.1  SCALAR  CASE 

Wbea  both  the  parameter  to  be  estimated  aad  the  measnrement  are  scalars,  the 
CRLB  reduces  to: 


EoYiroaraMt  li 


(D-15) 


where  p‘y(n)  deaotes  the  deriTative  of  P/v(<* )  respect  to  n  ,  and  /'(e)  denotes  the 
deriTatire  of  /  (e )  with  respect  to  e . 


EnvironniMt  St 


(D-IO) 


where  p'(e  )  denotes  the  deriTatire  of  p(e )  with  respect  to  e .  Note  that  the  fuDcthnal 
form  of  /  inflttences  only  the  Tirst  term  in  the  brackets,  whereas  the  a  priori  d>tribu- 
tioD  influences  both  terms. 


D.8.a  VECTOR  CASE 


One  may  obtain  similar  but  more  complicated  expressions  for  the  CRLB  for  the 
rector  case.  We  note  that  the  functional  form  of  /  influences  J  only,  whereas  the  a 
priori  distribution  influences  both  J  and  K . 
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D.94  nD  GAUSSIAN  CASE 


Here  ^  ■■  (Ni,  .  .  . ,  Nm)  n  t  vector  of  indepeadeBt  oad  ideoticillj 
(IID)  GouMiaa  r»Bdoa  variBblet,  tad  (ia  eaTiroBBeat  2)  A  i>  • 

IID  GoiufiaB,  ).  Tlie  CRLB  becomca: 


EaviroamMit  It 


Eoviroamrat  St 

£M. -•«.fl>|7  +  'fU'  . 

/.U)-5^  /.U)| 

K.f  —  t,,  .  ((,,  -  0  ,  I  y  »ad  —  1) 

(A  "  — /  ,  /  —  identity  matrix) 

<^4 


la  the  scalar  case  with  L  «  A/  »  1,  these  reduce  to 


distributed 

•  •  ,  )  is 

(D-17) 

(D-18a) 

(D-I8b) 

(D-IQ) 

(D-20) 

(D-213) 

(D-21b) 
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EnvlroamMii  It 


£0^  -.C|> 


Enviroamsot  St 


CM  > 


ev(a  )n/»^  +  i/»i 


or 


E\A  -A\*  ^  _ l_ 


<ri  “  Ey  (A  +  \ 


D.a.4  UNEAR  HD  GAUSSIAN  CASE 


Here 


fi  -  fa  +  iy  , 

where  F  a  M  X  L  mstrix;  A  ^<1  ^  ve  colamo  rectors;  sod 
iodependent  sod  HD  Goossioa  m  before.  In  this  case,  the  CRLB  becomes: 


Eariroamaat  It 


"  -2  L  fk  •‘hi  "  - -j - 


Eavlroamaat  St 


E{{Ai-A,n>[7  ^K]-: 


Fj_F 


H 


(D-22) 

(D-23) 

(D-24) 

(D-25) 
,  iV  are 

(D-26) 

(D-27) 

(D-28) 
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y  -  / 


(D-a9) 


F^F 

9n 


—  — 


(D-30) 


Additionally,  it  is  known  tkat  thcrs  exist  linear  estimators  for  which  the  bounds  are 
tight,  i.e.,  for  which  equality  holds. 


If  F  is  an  orthogonal  matrix,  i.e,  its  rows  are  orthonormal,  its  columns  are  ortho* 
normal,  and  F'^  ~  F  ^ ,  then 


Environment  It 

(D-31) 

Environment  St 

If  F  13  orthot^onal  except  for  a  scale  factor,  i.e.,  its  rows  are  orthogonal,  with  norm 
c  ,  then  f  *  F  I ,  and 


Environment  It 


(D-33) 


Environment  St 


E{{Ai  -A,)»)> 


l  +  c»(ri/w^ 


(D.34) 


If  F  is  0.  with  Fa  »  ,  as  it  wonid  he  if  it  represented  a  simple  weighting 

of  the  compopents  of  ^ ,  then 
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Envlroomrat  It 


<fN 


(D-3S) 


Eovironmant  St 


E{(Ai-Aif\> 


1+ 


(D-36) 


D.3  CRAMER.RAO  BOUNDS  IN  THE  PRESENCE  OP  AMBIGUITY 


Consider  the  additiTc  noise  estimstioi  problem  where 

+  ^  .  (D-37) 

If  /  is  not  one-to^ne,  we  ssy  thst  there  is  smlifwlf  in  the  function  /  sad  the  iness- 
urement  H  .  In  such  casee,  one  esanot  expect  any  estimator  to  perform  well.  Unfor* 
tunely,  howerer,  the  Cmmef>Rao  bound  may  not  reflect  this,  i.e,  it  may  giTe  a  rery  low 
MSE  eren  though  the  actual  MSE  b  rather  laige.  We  illustrate  thb  phenomenon  in  the 
scalar  case  and  show  that  it  b  not  a  problem  in  the  linear  Gaussian  case. 

Example  li  Scalar  Parameters 

One  may  see  from  Eqj  (D-15)  and  (D*I0)  and,  for  the  Gaussian  case,  Eq.s  (D-22) 
and  (D<23)  that  the  CRLB  will  be  the  same  for  any  function  /  whose  derivative  has  the 
same  magnitude  as  that  of  /  .  To  get  a  ckarer  pbture,  consider  the  situation  where  / 
has  a  continuous  second  derivative  but  b  not  one>U>oae.  Then  /  has  intervals  where  it 
b  constant  (/ '  »  0)  and/or  pairs  of  intervab  ,  /.  such  that  /  increases  on 
(/  '  >  0)  and  /  decreases  through  the  same  range  on  /.(/'<  0).  The  CRLB  will  be 
made  large  (and  appropriately  so)  by  the  intervab  where  /  '  »  0.  On  the  other  hand,  it 
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will  not  be  affected  by  the  exiitence  of  pain  of  interrab  where  /  increases  and 
decreases,  respectirely,  through  the  same  range,  i.e.,  it  will  be  as  small  as  the  CRLB  for 
the  nonambigttons  function 

/(•)-  /  .  (D-38) 

-oo 

for  which  one  expects  there  are  eatimaton  with  much  lower  MSE  than  for  /  .  Thus,  we 
conclude  that  in  the  presence  of  ambiguity,  the  CRLB  can  only  be  expected  to  gire  a 
reasonable  lower  bound  to  the  MSE  for  the  estimation  problem  inTolring  /  . 

As  a  concrete  example,  compare  the  estimation  of  A  based  on  either  but  not  both 
of  the  following  two  measurements 


,  (D-39) 

where  /  i(e )  «  s^  /  {(s )  «  |  e’  | ,  A  and  N  are  Gaussian,  independent  and  ij(0,  9^ ), 
f}(0,  9jii),  respectirely.  There  is  obrious  ambiguity  in  the  Rj  measurement.  For  either 
measurement,  the  CRLB  gires 


E[{A  -A)»l  > 


(D-^0) 


This  is  a  reasonable  bound  for  estimation  based  on  Ry  (It  tends  to  0  as  (t^v  ->  0;  it 
tends  to  oj  as  -•  oo.)  It  is  not  reasonable  for  estimation  based  on  ^2-  Indeed  f 
(Tjv  »  0,  then  the  best  estimate  for  A  based  on  R^i*  A  ^  E[A  \  /?2l  *  0  with  MSE 


(T ? ,  whereas  the  CRLB  lower  bound  is  aero. 


Finally,  consider  estimation  based  on  R^,  when  the  density  of  .t  is  replaced  by 
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ExunpU  3t  Linev  IID  GaoMian  Caae 
Here, 

-  FA  +  ^  ,  (D.42) 

where  fi  —  (/?„  .  .  .  ,  Ruf,  A  “■  (^  t,  •  •  • » )^,  —  (/Vi,  •  •  ,  ,  F  is  an 

M  X  L  matrix,  and  A  *Bd  ^  are  IID  GaoMiao  i|(0,  erj)  and  i|(0,  (rjii),  respectively  and 
independent  of  each  other.  For  ^lIMSE  estimator,  the  CRLB  bound  (10)  is  tight  and 
gives 

r  r  1 

-^+-1  (D-431 

L  »»  J. 

If  F  has  rank  less  than  L  (for  example,  if  M  <  L),  there  will  be  ambiguity  in  ^ . 
Nevertheless,  as  mentioned  in  Section  D.2.4,  the  CRLB  is  tight. 

As  a  concrete  example,  suppose 

F  —  A,  +  Aj  +  N  ,  (D-4-1) 

so  that  L  —  2,  iW  —  1,  F  [I  l|.  Then 
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1  +  <^4 

^N  ^A  ^N 

and  for  the  MMSE  estimator  it  is  well  known  that 


E[{A^-A,i^~E[{A,-A,f] 


fi(£N_+_^) 

ojlt  +  lai 


Note  that  as  -»  0,  MSE  -*  ell  2,  and  as  -»  oo,  MSE  -•  w j . 


(D.45) 


(D-40) 


(D-47) 
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,  Appendix  E 
^*5  (g(x))FOR  COMPLEX  OBJECTS 

In  this  appendix  we  generalize  the  expression  for  the  gradient  of 
the  summed  objective  function  to  Include  objects  and  object  estimates 
that  are  complex  valued. 

Recall  that  the  sunmed  objective  function  Is  defined  as  the  sum  of 
a  generalized  object-domain  error  metric  and  the  Fourler-domaln  error 
metric: 


e 


2 

s 


where 


(E-1  ) 


(E-2  ) 


and 


[  l6(u)|  -  lF{u)l]^  (E.3  ) 
u 

Notice  the  summed  objective  function  Is  Implicitly  a  function  of  the 
real  and  Imaginary  parts  of  the  pixel  values  for  the  latest  estimate, 
g(x).  We  therefore  treat  the  real  and  Imaginary  parts  of  each  pixel  as 
distinct  parameters  that  can  be  adjusted  in  order  to  minimize  e^  .  We 
express  the  real  and  iMgInary  parts  of  the  latest  estimate  as 
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g(x)  -  a(x)  *  ib(x). 


(E-4) 


We  define  the  gradients  to  be: 


7«|(g(«))  ■  E 


aa(xj) 


R  ^ 

J  3b(Xj 


(E-5) 


R  I 

where  and  Vx  are  orthogonal  unit  vectors  associated  with  the  real 
J  J  th 

and  Imaginary  parts  of  the  J  pixel.  The  first  partial  derivative  In 
Eq.  (E-5)  may  be  separated  Into  two  terms: 


3a(xj)  3a(xj)  3a(xp 

For  the  moment  we  examine  the  first  term  In  (E-6): 


u 

•  2N-2  ^  (|G(u)|  -  |F(u)|) 


It  Is  easy  to  see  that 


2Ti:rt 


br[' 


G(u)  e 


IZttu-x  ./N 


>  C 


‘1 


(E-6) 


(E-7) 

(E-8) 


(E-9) 
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Substituting  Eq.  (E-9)  into  Eq.  (E-8)  yields 


N-2  ^  (|G(u)|  -  |F(u)|)  MhIA 


i27iu*)c  /N 


TgMT 


c.c.) 


«  (g(Xj)  -  g'(xj))  +  (g*(xj)  -  g**(Xj)) 


*  2(a(Xj)  -  a'(Xj)). 


(E-IO) 

(E-11) 

(E-12) 


This  result  is  consistent  with  the  result  quoted  for  real-valued 
objects.  A  parallel  derivation  gives 


(E-t3) 


We  now  return  to  the  second  term  in  Eq.  (E-6): 

Y,  l9(x)l^ 


.  3 

3a(x J  9a(x.) 


xcS' 


Similarly, 


^  xcS' 

2a(x.)  .  X.  c  S- 

0  ,  X  .  c  S 

J 


3e  2  f2b(x.)  .  X.  e  S’ 

0  _  J  J 

3b(x .)  " 

0  .  X  e  S 


(E-14) 


(E-15) 
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Let 


S'(x) 


1.  X  e  S' 

(  0.  X  e  S 


(E-16) 


Now  Eqs.  (E-14)  and  (E-15)  may  be  expressed  more  conveniently: 


5IT^ 


2a(x.)S'{x.) 


(E-17) 


and 


=  2b(x.)S'(x.) 


(E-18) 


Collecting  these  results,  we  have; 


»  2[a(xj)  -  a'(xj)]  2a(xj)S'(x.) 


(E-19) 


t3(x*)  “  2[b(Xj)  -  t>‘(Xj)]  +  2b(Xj)S'(Xj) 


(E-20) 


Equations  (E-19)  and  (E-20)  may  be  combined  to  form  a  complex  gradient 
Image. 
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Gradient  Image  =  2[g(x)  -  g'(x)]  ♦  2g(x)S'(x) 


(E-21) 


The  extension  to  complex-valued  objects  still  requires  only  2  FFTs  to 
compute  the  gradient. 
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,  Appendix  F 
Ve^  (g(x))  FOR  REAL  OBJECTS 

Recall  that  the  object-domain  error  metric  Is  Implicitly  a  function 
of  the  Input  estimate  g(x)  and  Is  given  by: 

•,^(3(x))-  2[9'(x)]^  (F-') 

Its  gradient  with  respect  to  the  input  pixel  values  may  be  written: 

X  WITT  "i  (F-2) 

J-l  ■’ 

Where  Vj  Is  a  unit  vector  In  the  direction  of  the  parameter  g(xj)  In 
parameter  space.  We  focus  now  on  the  partial  derivative  that  appears  in 
Eq.  (F-2): 


Substitution  of  the  Fourler-domain  expression  for  g'(x)  gives: 


/ 


g'(x) 


35nj) 


^  G'(u) 


Recall  that 


(F-4) 


G'(u) 


(F-5) 


176 


giving 


xeS' 


gi2Tr'j-x/N 

4 


z 

xeS' 


g‘(x) 


^  [F(u)|e^2nu-x/N 


(F-6) 


In  order  to  evaluate  the  partial  derivative  In  Eq,  (F-6)  we  need 

m- 


expressions  ear  end 

SsTx-l 


”  39TfjT  Z 


-12nu-x/N 


=  e 


•12Tru*x/N 


i2Ttu-Xj/N 


The  dcrivetlon  for  repo  ires  the  result  In  £q.  (F-7); 

G(u) 


9  i G( u )  [  ,  _ 1  ^ _ 


3qix.i  rZUIJT 


JTTOTT 


G*(u) 


-i2Tiu-x./N 

G*(u)  e  J  +  C.C 


■] 


(F-7) 


(F-8) 
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where  C.C.  stands  for  complex  conjugate  of  the  explicit  term.  Using  the 
results  In  (F-7)  and  (F-8)  and  with  some  algebraic  manipulation  the 
partial  derivative  In  (F-6)  becomes 


a  G(u) 
3g(Xj)  |6(u)| 


-12iru*x./N 

G*(u)  e  J  -  C.C. 
2|6(u}I  6*(u) - 


(F-9) 


Substituting  Eq.  (F-9)  back  Into  Eq.  (F-6)  yields 
2 


J  xcS'  li 


-12itu*x./N 

G*(u.)  e  J  -  C.C. 

|fi(u)I6*(u) - 


i2iTU'x/N 


(F-10) 


By  changing  the  order  of  summation  we  have 

Y,  9'(x) 

xeS'  (F-11) 

At  this  point  It  Is  convenient  to  define  the  characteristic  function  of 
the  complement  of  the  support  S  as  follows 


3e 


F(u) 


-i2iru-x./N 

^  -  C.C. 


IG(u)lG*(u) 


S'(x) 


1  .  X  c  S’ 
0  •  X  c  S 


(F-12) 


The  second  suMMtIon  in  (F-11)  may  now  be  rewritten 


y  g’(x) 
xeS* 


i2Tru’X/N 


5]  S’(x)g’(x)  e 


i2nu- x/N 


(F-13) 


The  error  In  the  output,  g^(x),  consists  of  that  component  of  the  output 
that  violates  the  support  constraint: 


g,(x)  ■  S'(x)g’(x) 


(F-14) 
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Th«  suMMtIon  in  Eq.  (F-13)  has  the  form  of  a  forward  DFT: 

^  9,(x)  ^  5^(x) 

X  X 


■  Gg(-u) 


(F-15) 


whtre  G^(u)  Is  th«  OFT  of  g^(x).  Tht  total  partial  dtrivativa  may  now 
ba  wrlttan: 


2 

'*0 

3g(Xj) 


|F(u)|Gg{-u)  -12itu-x./N 

|6(u)|i*(u)  “*<“>  * 

■I 


.2  r'  l^(u)|G  (-u)  -IZttu-x./N  j  _  lF{u)lG,(-u)G(u)  i2~u  x  /N 

■ "  L  ii(.)i  '  '  - "  Z  rroiiMu) — '  ' 

“  “  (F-U) 

Using  Eq.  (F-5)  In  tht  stcond  term: 


3e2  .  _  |F(u)|G  (-u)  -i2T»u-x./N 

Tgfl-T-"  Z  -iSMl  * 


G'(u)G  (-u)  iZ^u-x./N 

S*TuT  " 

•  (F.17) 


Remarkably  both  of  these  terms  have  the  general  form  of  OFTs.  In  order 
to  combine  both  of  these  terms  Into  a  single  Inverse  DFT  we  perform  a 
sign  change  of  variable  on  the  Fourier  vector  In  the  first  term.  The 
net  result  Is: 


K 


x-2 


^  f|F(-u)|G^(u) 

2-  iG(-u)  1 

u  t 


G'(u)Gg(-u)  IZiru-x./N 


(F-18) 


179 


Finally,  by  appealing  to  the  Hermltlan  property  for  Fourier  transforms 
of  real  functions  we  may  make  the  following  substitutions: 


|F(-a)|  «  |F(u)| 

|G(-u)|  «  |G(u)| 

G^(-u)  »  Gg*(u)  (F.19) 


to  get  the  final  result: 


3«o  -2  r* 

39f5TT"'  I  |i(.)l 

J  U 


G‘(u)G^*(u) 

G*(u} 


12nu*x./N 

e  ^ 


(F-20) 
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.  Appendix  G 

^0  (g(x))  for  complex  objects 

The  derivation  of  the  gradient  of  the  (g(x))  objective  function 
for  the  case  of  complex  objects  closely  follows  that  for  real  objects 
(Appendix  F).  Therefore,  we  Just  highlight  this  derivation,  calling 
attention  to  significant  differences.  Ue  begin  by  acknowledging  that 
the  complex  case  admits  twice  as  many  parameters;  namely  the  real  and 
Imaginary  parts  of  each  Input  pixel.  We  denote  the  real  and  Imaginary 
parts  of  the  Input  function  as  follows: 


g(x)  ■  a(x)  ♦  1  b(x). 


(G-l) 


We  define  the  gradient  to  be: 


VeoCg(x)) 


{G-2) 


R  I 

where  Vj  and  v^  are  orthogonal  unit  vectors  associated  with  the 
parameters  of  the  real  and  Imaginary  parts  of  the  pixel  at  location  Xj. 
The  partial  derivative  of  the  objective  function  with  respect  to  the 
real  part  of  a  pixel  value  may  be  written 

ae2 

J  J  xeS' 


C.C. 

(G-3) 
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The  extra  complex-conjugate  term  appears  because  the  output  function, 
g'(x),  can  assume  complex  values.  The  partial  derivative  with  respect 
to  the  Imaginary  part  is  similarly  found: 

2 

*  C.C. 
(G-4) 

With  simple  algebraic  manipulations  the  following  useful  Identities  may 
be  verified: 


3e 


-2 


I 

xeS' 


g'*(x) 


,12itu*x/N 


3 

3R? 


G(u) 


{G-5) 

(G-6) 

(G-?) 

(G-8) 


With  the  aid  of  Eqs.  (G-5)  thr*u  (G-8)  we  deduce 


-iZTTux  ./N 

e _ ^ 

2|G(u)1G*(u) 


C.C. 


(G-9) 
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and 


G 


•12ffU-x./N 


^|G(u)iG*(u) 


-  C.C 


(G-10) 


When  Eq.  (G-9)  Is  substituted  Into  Eq.  (G-3)  and  the  order  of  summation 
Is  Interchanged  we  get: 

2 


3a (Xj)  “ 


N-2 

T 


-iZffu-Xj/N 


^  S-(x)g-*(x) 


x/N 


C.C. 


(G-  1 1 ) 

were  S'(x)  Is  the  characteristic  function  of  the  set  S',  as  before. 
Recall  the  object  and  Fourler-domain  expressions  for  the  error  image: 


g^(x)  -  S'{x)g'(x) 

Ge(u)  -  Y,  S'(x)g'(x)e"^^’'“**^'' 


(G-12) 

(G-13) 


We  may  therefore  write 


Ge*(u)  -  Y  S'(x)g'*(x)e^^’'“’*^'’ 


(G-14) 


which  may  be  substituted  Into  (G-11)  to  get: 


3a(Xj) 


■  i  -i2TTU-x,/N 

^-2^  |F(u)|Gg*(u)  jG*(u)  e  J 

'~r  h 


16(u)|G*(u) 


C.C. I 


+  C.C.  (g-  15) 
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By  similar  steps  we  find 


n-2 

T 


F(u)|G  * 

_ e 


(u)|lG*(u)  e 
— 


-12iTU-Xj/N 


:.c| 


+  C.C.  (G-16) 


Explicitly  writing  out  the  complex-conjugate  terms  In  Eqs.  (G-15)  and 
(G-16)  and  performing  a  few  additional  manipulations  we  produce  the 
final  result. 


^^0 

3b(Xj) 


Re 


""I 

u 

t 

■|F(-u)|G^*(-u) 

G'(u)Gj*(u)] 

|G(-u)| 

• 

~5*(u)  J 

u 

’F(-u)Gg*(-u) 

G'(u)G^*(u)'| 

|6(-u)|  * 

G*(u)  ^1 

i2TTu-Xj/N 


(G-17) 


i2iTu*x  ./N 

J 


(G-18) 


It  Is  gratifying  that  Eq.  (G-17)  Is  consistent  with  Eq.  (F-18)  which  Is 
the  equivalent  partial  derivative  for  real  objects  only.  It  Is  worth 
mentioning  that  because  the  summation  arguments  In  Eqs.  (G-17)  and 
(G-18)  differ,  an  additional  FFT  Is  required  In  the  computation  of  the 
gradient  for  complex  objects.  The  total  number  of  FFTs  (forward  or 
Inverse)  needed  Is  Increased  to  five. 
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